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FOREWORD 


The  goal  of  Range  Safety  is  the  prevention  of  injury  to 
personnel  or  damage  to  property  by  taking  all  reasonable 
precautions  consistent  with  operational  requirements.  This 
is  dependent  not  only  on  precautions  taken  in  the  preparation 
of  a missile  or  vehicle  launch,  but  in  the  ability  of  the 
Range  Safety  Officer  (RSO)  to  maintain  surveillance  during 
flight  to  insure  compliance  with  established  safety  criteria. 
To  maintain  this  necessary  surveillance,  the  RSO  must  have 
at  his  disposal  information  depicting  performance  of  the 
missile  and  possible  impact  locations  for  comparison  against 
predetermined  destruct  criteria  as  well  as  assurance  that 
his  e safety  system  is  in  operational  condition  at  all 

t renting  this  information  in  a manner  that  allows 

ud  quick  understanding  of  significant  data  will  vary 
due  to  missile  dynamics  and  test  range  geometry.  Therefore, 
there  is  no  one-best  filter  for  all  applications. 

A survey  of  seven  of  the  test  ranges  represented  in  the 
Range  Safety  Group  (RSG)  of  the  Range  Commanders  Council 
(RCC)  was  made  to  determine  the  type  of  filter  systems 

currently  in  use  at  the  various  ranges.  This  information  is 

. •£ 

presented  in  this  document  as  an  aid  to  all  ranges  in  deter- 
mining which  systems  may  have  mejfit  for  their  application 
and  to  provide  some  insight  into,  future  applications  of 
filter  systems. 
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1.0  INTRODUCTION 


1.1  Flight  Safety  Data  Working  Group 

1.1.1  The  formation  of  the  Flight  Safety  Data  Working 
Group  (FSDWG)  was  originally  considered  at  the  16th  meeting 
of  the  Inter-Range  Missile  Flight  Safety  Group  (IRMFSG). 
Member  ranges  were  requested  to  select  nominees  for  this 
working  group..  At  the  17th  meeting  of  the  IRMFSG,  held  28- 
29.  April  1964,  the  FSDWG  was  established.  The  purpose  of 
this  working  group  was  to  assist  the  IRMFSG  by  studying, 
making  recommendations  to  the  IRMFSG  chairman,  and  promoting 
the  exchange  of  information  between  member  ranges  concerning 
generating,  processing,  and  displaying  flight  safety  data 
for  the  Flight  Safety  Officer  (FSO). 

The  second  and  subsequent  meetings  of  the  FSDWG 
were  held  at  the  following  ranges:  Western  Test  Range, 
October  1964;  NASA  Houston,  June  1965;  Colorado  Springs, 

March  1965;  Eastern  Test  Range,  June  1966;  White  Sands 
Missile  Range,  November  1966;  Pacific  Missile  Range,  May 
1967;  Eastern  Test  Range,  May.  1968;  Wallops  Station,  October 
1968;  Air  Force  Satellite  Control  Facility,  April  1969;  Air 
Force  Electronics  Systems  Division,  October  1969. 

1.1.2  From  April  1964  through  June  1966,  the  prime 
discussion  of  the  FSDWG  centered  around  the  development  of 
the  various  impact  prediction  systems  being  used  at  the  RCC 
member  ranges.  At  the  November  1966  meeting  (6th  meeting), 
the  group  began  a detailed  investigation  into  the  filtering 
techniques  employed  at  these  ranges.  A detailed  investigation 
into  the  Kalman  technique  was  included.  At  the  9th  meeting 

at  Wallops  Station  in  October  1968,  it  was  agreed  that  each 
range  would  develop  a subroutine,  for  exchange  with  other 
ranges,  to  be  utilized  for  comparison  of  filter  techniques. 
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This  exchange  occurred  at  the  10th  meeting  of  the  FSDWG  and 
emphasis  was  then  placed  upon  the  development  of  a report  of 
"Status  of  Range  Safety  Filter  Systems,"  Task  6-25,  issued 
at  the  25th  IRMFSG  meeting  of  the  IRMFSG,  25-26  September 
1968. 

1.1.3  An  ad  hoc  committee  was  established  to  revise  the 

original  document  and  to  add  an  ADTC  section  in  March  1976. 

1.2  Report 


1.2.1  Task  6-25  was  given  to  FSDWG  in  an  effort  to  make 

a survey  of  real-time  filtering  and  editing  techniques 
presently  in  existence  and  publish  the  results  in  an  IRMFSG 
document.  This  survey  was  to  include  system  characteristics, 
advantages,  disadvantages,  etc. 


1.2.2  At  the  FSDWG  meeting  held  in  October  1969,  at  the 
Air  Force  Electronic  Systems  Division,  a proposed  outline  of 
the  task  was  developed.  This  outline  was  arranged  to  provide 
four  specific  chapters  covering  the  following  areas:  (1)  a 
brief  introduction  into  the  background  of  the  FSDWG  and  Task 
6-25  (including  the  purpose  of  the  report  and  the  method  of 
preparation);  (2)  a discussion  of  filtering  as  applied  to 
flight  safety,  the  instrumentation  systems  being  utilized, 
and  the  computation  and  display  systems  available  at  the 
various  ranges;  (3)  a brief  and  general  outline  of  the  types 
of  filters  being  utilized  with  some  descriptive  information 
concerning  their  specific  characteristics,  and  (4)  the 
specific  filters  in  use  at  the  ranges  surveyed,  how  and  why 
they  were  developed,  the  mathematics  and  operation  of  the 
filters  and  their  outstanding  features  and  limitations. 

1.2.3  Initial  discussions  for  accomplishing  this  task 
indicated  a need  for  general  background  information  on 
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filtering  systems  and  detailed  information  on  specific 
filters  in  use  at  the  RCC  ranges.  Each  range  was  requested 
to  supply  detailed  information  on  filter  systems  in  use,  why 
the  system  was  developed,  how  it  is  used,  theory  of  the 
system,  outstanding  features  and  limitations,  and  types  of 
instrumentation  data  filtered. 

2.0  DISCUSSION 

2.1  The  purpose  of  filtering  data,  for  range  safety  systems, 
is  to  remove  unwanted  noise  from  the  measured  data  in  order 
to  determine,  as  accurately  as  possible,  the  space  position, 
velocity  and  acceleration  of  the  target.  There  are,  no 
doubt,  as  many  filtering  techniques  being  used  as  there  are 
persons  evaluating  the  results.  However,  with  the  proper 
computations,  the  results  provide  the  FSO  with  information 
which  can  be  readily  used  in  determination  of  the  actual 
trajectory.  In  order  to  accomplish  this  objective,  data 
measurements  from  such  instrumentation  as  radars  are  entered 
into  a computer  system,  normally  at  high  input  rates  such  as 
10-20  samples  per  second.  These  data  must  then  be  differen- 
tiated to  obtain  velocity  and.  acceleration  for  use  in  calcula- 
tion of  predicted  impact  point  and  output  to  displays  for 
use  by  the  FSO.  However,  since  all  instrumentation  systems 
produce  erroneous  data,  errors  must  be  removed  before  differ- 
entiation and  prediction  can  occur-.  The  computer  system  is 
also  limited  in  speed  and  memory,  therefore  limiting  the 
ability  to  statistically  examine  the  data  samples  in  real  time 
These,  in  general,  are  the  problems  associated  with  data 
filtering  for  flight  safety.  The  methods  used  at  the  various 
ranges  surveyed  are  described  more  fully  in  Section  4.0  of 
this  document. 

2.2  Instrumentation  Systems 

2.2.1  A wide  variety  of  instrumentation  systems  is  used 
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to  provide  the  FSO  with  information  required  to  perform  the 
decision-making  process.  This  section  is  concerned  with 
those  instrumentation  systems  which  provide  data  for  input 
into  the  real-time  computer.  Historically,  pulse  radars 
have  been  the  more  popular  of  the  real-time  range  safety 
instrumentation  systems.  These  types  of  radars  primarily 
measure  the  position  vector  only,  i.e.,  range  (R),  azimuth 
(A)  and  elevation  (E).  Recently,  doppler  rate  measurement 
devices  have  been  added  to  some  of  these  pulse  radars  to 
allow  direct  measurement  of  range  rate.  Most  of  these 
devices,  with  the  exception  of  the  FPS-16  rate  kits,  can 
operate  with  or  without  a coherent  beacon;  however,  more 
accurate  range  rate  measurement  is  obtained  with  a coherent 
beacon.  Coherent  beacons  are  in  the  development  stage  and 
their  usage  is  expected  to  increase  in  the  future. 

Continuous  wave  (CW)  systems  have  found  usage 
where  extremely  accurate  position  and  velocity  measurements 
are  required.  These  can  serve  the  dual  purposes  of  guidance 
system  evaluation  and  range  safety.  For  example,  GERTS 
provides  a direct  measurement  of  three  components  of  the 
velocity  vector  (R,  P,  and  Q)\  These  systems  are  expensive 
to  purchase  and  operate,  and  normally  require  expensive, 
special  purpose  transponders  onboard  the  vehicle.  There  are 
other  varied  types  of  measurement  systems  in  use  at  the 
ranges  surveyed.  Some  of  these  measure  position  through 
interferometer  techniques;  others  measure  range  rate  from 
doppler  frequency  shifts.  Many  of  these  systems  are  mobile 
sensors.  Several,  in  combination,  may  be  required  to  produce 
either  a position  vector  for  differentiation  or  a velocity 
vector  for  combination  with  position  measured  from  other 
sources . 


Data  available 
inertial  and  radio,  have 
computations.  Most  radio 


from  missile  guidance  systems,  both 
gained  in  usage  in  range  safety 

guidance  systems  are  either  pulse 
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radars  or  CW  systems  as  discussed  above.  Data  from  inertial 
guidance  systems  have  found  use  primarily  due  to  the  lack  of 
sufficient  accuracy  obtained  through  either  pulse  or  CW 


systems.  The  use  of  telemetered  inertia]  guidance  data  for 
range  safety  purposes  is  questionable  from  the  standpoint 
that  the  system  which  is  guiding  a malfunctioning  missile  is 
also  being  used  in  the  computations  to  determine  whether  it 
should  be  destroyed.  Experience  has  shown  that  the  earliest 
indication  of  a malfunction  is  often  obtained  through  on- 
board telemetered  measurements. 

2. 2. 1.1  Instrumentation  systems  that  are  currently  installed 
and  operating  at  the  ranges  surveyed  are  identified  in  Table 
1. 

2. 2. 1.2  Few  new  instrumentation  systems  that  can  be  used 
for  range  safety  purposes  are  planned  in  the  immediate 
future.  For  the  most  part,  existing  systems  will  continue 
in  use  and  will  be  representative  of  future  systems.  Efforts 
in  radar  calibration  are  expected  to  result  in  improved 
tracking  accuracies.  In  addition  to  inertial  guidance  data, 
increased  range  safety  use  of • telemetered  parameters  such  as 
guidance  commands,  nozzle  deflection  angles,  and  pitch,  yaw, 
and  roll  angles  and  rates  can  be  expected.  Data  from  external 
tracking  systems  will  be  combined  in  the  real-time  computer 
with  telemetered  guidance  data  and  missile  body  parameters 
to  provide  additional  data  for  range-safety  decision  making. 
With  upgraded  computers,  analytical  models,  and  display 
systems,  range  safety  visibility  into  real-time  missile 
performance  should  improve  correspondingly. 

To  provide  accurate  post-flight  trajectory  data 
for  Trident  missile  launches  on  the  Eastern  and  Pacific 
missile  ranges,  the  missiles  will  be  equipped  with  special 
translators  that  operate  in  conjunction  with  orbiting  satel- 
lites and  ground  stations  referred  to  as  "pseudo  satellites." 
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TABLE  I 


Installed  Range  Instrumentation  Systems 


DATA  TRACKING 

DESCRIPTION  OBSERVED  MODE*  MANUFACTURER  RANGE** 

A.  Pulse  Radars 


MOD  II 

R,  A,E  . 

T , S 

REEVES 

AN/FPS-16 

R , ,A , E , ( R ) 2 ’ 3 

T,  S 

RCA 

AN/FPQ-6 

R , A , E , R 

T , S 

RCA 

AN/TPQ- 18 

• R , A , E , R 

T , S 

RCA 

AN/MPS-25 

R , A , E • . 

T , S 

, RCA 

MPS-19 

R , A , E . 

T , S 

REEVES 

SPANDAR 

R , A , E , R 

T , S 

AN/FPQ-10 

R , A , E 

T , S 

Sperry 

AN/FPS-13 

R , A , E . 

T , S 

RCA 

MPS-36 

R , A , E , R 

T , S 

RCA 

AN/FPQ-14 

R , A , E*** 

T , S 

RCA 

ALCOR 

R,A,E**** 

T , S 

RCA 

B.  CW  Systems 

GERTS 

R , A , E , R , P ,Q 

T 

GE 

X,Y,Z 

T 

C.  Guidance  Systems 

Inertial  Position  & Telemetry  1,2,4 

Velocity 
Components 


Radio 


R , A , E , 


2,6 


*T  = Transponder  or  beacon;  S = Skin  tracking. 

**1  = DET  1/SAMTEC , 2 = SAMTEC , 3 = WSMR,  4 = PMTC,  5 = Wallops, 
6 = KMR , 7 = ADTC . 

( )x  Denotes  additional  capability  existing  at  X range 
over  and  a we  the  standard. 

***  Data  observed  in  R,A,E,  but  output  in  geocentric, 
earth  fixed  EFGE#G  position  and  velocity. 

****ARPA  Lincoln  C-band  observables  radar. 
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By  adding  other  remote  ranging  stations,  the  system's  appli- 
cation will  be  extended  to  include  real-time  range  safety 
use.  Four  remote  stations  can  be  used  simultaneously  in  a 
time-division-multiple-access  mode  with  a master  ranging 
station.  The  master  and  each  remote  station  transmit  signals 
in  sequence  to  the  missile  where  the  signals  are  shifted  in 
frequency  by  the  missile  translator  and  retransmitted  to  the 
master  station.  Here  the  signals  are  processed  to  extract  a 
range  sum  and  range-rate  sum  for  each  station.  Missile 
position  and  velocity  components  are  then  computed  from  these 
sums . 

2.2.2  Range  safety  computation  and  display  equipment, 

both  that  in  use  and  that  planned  for  installation  in  the 
future,  represents  a variety  of  well  known  companies. 
Supervisory  software  techniques  employed  with  computational 
equipment,  while  quite  varied  among  the  ranges  surveyed, 
generally  fall  within  two  categories:  (1)  dedication  of  a 
computer  system  to  the  solution  of  range  safety  related 
problems  is  noted  in  the  next  two  subparagraphs  as  dedicated 
(D);  and  (2)  operation  within  a multiprogrammed  and  possibly 
multiprocessor  environment  which  allows  other  unrelated  data 
processing  to  continue  concurrently  with  the  real-time  range 
safety  calculation  is  noted  in  the  next  two  subparagraphs 
as  shared  (S). 

Display  equipment  in  use  for  presenting  range 
safety  calculations  are:  plotters  (P),  digital  displays 
(DD),  cathode  ray  tubes  (CRT),  teletypes  (TTY),  and  lighted 
matrices  (LM).  Lighted  matrices  generally  are  used  to 
present  instrumentation  and  measurement  data  status. 
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2.2.2. 1 Presently  Installed  Systems 


Range 

Computer 

Computer  Usage 

Display  Equip. 

SAMTEC 

IBM-7044 

D 

P , LM 

DET  1 /SAMTEC 

CDC-3600/CDC  3100 

D 

CRT , LM , DD 

NASA/Wallops 

HW-625 

D , S 

DD , P , CRT , TTY , 

LM 

WSMR 

Univac  1108 

D , S 

P , DD , LM , TTY 

PMTC 

Univac  1230/1212 

D 

P , LM 

KMR 

CDC  7600/PDP  11-05 

S 

CRT , LM 

ADTC 

CDC-6600/ IBM  360-65/ 

D , S 

CRT , P 

PDP-15/PDP-11 


2. 2. 2. 2 Planned  Future  Systems 


Range 

Computer 

Computer  Usage 

Display  Equip. 

SAMTEC 

IBM  360-65  to  be 

replaced  by  IBM 

370-165 

D 

P , DD , CRT , LM 

DET  1 /SAMTEC 

CDC-CYBER 

74 

D , S 

CRT , LM , DD 

NASA/Wallops 

None 

P , TTY , LM , 

expanded 

use  of  CRT 

WSMR 

Upgraded  System 

Planned 

D 

P , DD , LM , CRT , 

TTY 

PMTC 

Upgraded  System 

Planned 

D 

CRT  ( Color ),P, 

DD , TTY , LM 

3.0  THEORY  OF  FILTERING 

3.1  The  basic  objective  of  filtering  is  to  remove  random 
noise  from  observations.  This  noise  may  disguise  valid  data 


to  a point  where  it  cannot  be  used  as  an  accurate  measurement. 
The  techniques  used  to  accomplish  this  may  vary  in  complexity 
from  the  averaging  of  data  samples  to  the  more  complex 
derivation  of  state  vectors  for  prediction  of  new  measurements 
and  recursive  updating  of  the  estimates.  The  most  common 
techniques  employed  in  the  past  have  been  the  development  of 
finite-length  least-square  filters.  With  the  improved 
computational  systems  the  filtering  of  data  may  be  accomplished 
by  more  complex  statistical  techniques. 

In  the  selection  of  a mathematical  technique  to  provide 
the  best  determination  of  position,  velocity,  and  acceleration 
one  must  always  consider  the  limitations  imposed  by  the 
computational  equipment.  This  means  that  the  technique  must 
be  adapted  to  the  solution  of  the  problem  as  well  as  the 
computing  equipment  to  be  used.  Such  adaptation  is  especially 
necessary  when  considering  a filtering  technique  for  use  in 
real  time  to  supply  inputs  to  range  safety  calculations. 

Restrictions  imposed  through  implementing  any  computation 
in  real  time  are:  (1)  computational  time  requirements  and 
(2)  computer  storage  requirements.  Obviously,  a sophisticated 
best  estimate  of  trajectory  scheme,  which  requires  64,000 
words  of  computer  storage  and  1 second  of  computational  time 
per  cycle,  cannot  be  implemented  on  a 32,000  word  storage 
computer  ten  times  per  second. 

3.2  Filtering  Techniques 

3.2.1  Data  filtering  or  smoothing  can  be  thought  of  as 

the  replacement  of  a curve  or  a sequence  of  points,  by 
another  that  is  in  some  sense  more  regular,  and  yet  whose 
ordinates,  for  any  abscissa,  are  changed  as  little  as  possible. 
The  irregularities  in  sequence  of  points  may  be  due  to 
errors  in  measurement.  If  theory  requires  the  theoretically 
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correct  points  to  lie  on  a given  curve,  one  may  apply  some 
method  of  curve  fitting,  possibly  least  squares.  If  not, 
one  may  select  arbitrarily  a simple  function,  possibly  a 
polynomial,  and  fit  it  by  least  squares.  If  the  purpose  is 
merely  to  obtain  a smooth  graph,  this  may  be  drawn  visually. 

A somewhat  more  sophisticated  method  is  to  take,  say,  5 
consecutive  points,  fit  a parabola,  and  replace  the  middle 
point  by  the  one  on  the  parabola.  The  next  parabola  requires 
four  of  these  points  and  one  new  one. 

3.2.2  Filtering  techniques  can  be  roughly  divided  into 

two  basic  groups  as  follows: 

(1)  Recursive  filters  may  operate  on  a stored 
array  of  previous  measurements  or  restrict  calculations  to 
the  latest  point.  The  principle  of  recursion  as  applied  to 
data  filtering  relates  the  previous  filter  output  to  the 
current  calculated  values.  Such  a recursive  filter  cannot 
effectively  operate  without  the  previously  calculated  values. 

This  implies  that  the  filter  must  be  supplied  an  approximation 
of  these  values  for  any  initial  filter  start.  Many  schemes 
for  supplying  these  initial  approximations  of  position, 
velocity,  and  acceleration  are  in  use  and  are  known  as 
filter  initialization. 

(2)  Non-recursive  filters  generally  operate  on  a 
stored  array  of  raw  data  points.  This  type  of  filter  requires 
only  that  the  array  be  filled  and  updated  as  each  new  point 

is  received.  No  history  of  previous  filter  output  is  required 
to  generate  the  current  values  of  position,  velocity  and 
acceleration . 

The  recursive  and  non-recursive  properties  of 
filters  are  sometimes  difficult  to  categorize  and  therefore 
examples  of  each  type  are  presented  in  the  following  subparagraphs. 
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3.2.2. 1 The  most  simple  non-recursive  technique  of  filtering 
(smoothing)  data  is  averaging  points.  This  can  be  a simple 
average  of  data  points  to  establish  a mean  or  midpoint  on  a 
curve  or  a more  complex  weighted  average  designed  to  establish 
the  latest  value  or  end  point.  Average  velocity  can  be 
determined  by  computing  the  average  of  position-first  differ- 
ences. Average  acceleration  can  be  determined  by  computing 
the  average  of  position-second  differences.  The  latest  or 
real-time  values  of  velocity  and  acceleration  can  be  approx- 
imated by  weighting  the  average.  Many  pre-filter  data  editing 
or  filter  initialization  schemes  use  averaging  techniques. 

3. 2. 2. 2 The  approach  of  filtering  (smoothing)  the  raw 
measurements  often  utilizes  a "classical"  least-squares 
technique.  This  technique  simply  utilizes  the  principle  of 
fitting  a given  number  of  data  points  to  a second-  or  third- 
degree  polynomial  so  as  to  minimize  the  sums  of  the  squares 
of  the  residuals  between  the  polynomial  and  the  measured 
data  points.  By  making  certain  assumptions  such  as  the 
number  of  data  points,  the  output  point,  constant  time 
intervals,  etc.,  the  least-squares  technique  can  be  resolved 
into  a finite  number  of  coefficients,  called  weights,  which 
can  then  be  used  to  scale  the  input  data  points  to  yield  the 
desired  output.  For  purposes  of  flight  safety  this  point 
must  be  related  to  the  most  recent  input  value  in  order  to 
avoid  costly  time  delays.  The  above  technique,  however, 
does  not  consider  such  important  information  as  vehicle 
capabilities,  gravity  models,  and  instrumentation  noise 
estimates.  In  order  to  introduce  these  effects,  the  above 
techniques  may  be  modified  to  make  use  of  this  information 
as  described  in  the  following  subparagraphs: 


3. 2. 2. 3  A variety  of  formulations  of  the  Kalman  filter 
exist.  They  all  provide  a means  for  the  optimum  use  of  a set 
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of  measurements  to  give  the  best  (i.e.,  least  squares) 
estimate  of  the.  state  of  a system  and  are  based  on  the  work 
of  Dr.  Rudolf  E.  Kalman.  The  Kalman  filter  may  be  expressed 
very  generally  as  a set  of  matrix  equations,  but  the  content 
of  the  matrices  will  vary  depending  on  the  system  to  which 
the  filter  is  applied.  Many  filters,  however,  are  erroneously 
referred  to  as-  Kalman  types  which  do  not  contain  the  essential 
optimal  data  weighting  feature. 

Those  parameters  of  the  system  of  interest  which 
adequately  describe  its  present  dynamic  conditions  are 
defined  as  the  state  variables  of  the  system.  For  example, 
the  state  variables  for  a missile  in  flight  might  be  the 
three  components  each  of  position,  velocity,  and  acceleration. 
These  nine  parameters  would  then  be  considered  the  coordinates 
of  a nine-dimensional  space,  and  the  condition  of  the  missile 
would  be  described  by  a nine-element  state  vector.  Those 
parameters  of  the  system  which  are  measured  are  considered 
to  be  the  elements  of  a measurement  vector.  For  example, 
with  a single  conventional  radar  there  would  be  a three- 
dimensional  measurement  vector  with  elements  of  range, 
azimuth  and  elevation. 

The  state  of  the  system  is  assumed  to  change  in 
some  predictable  manner  so  that  future  states  may  be  predicted 
from  the  current  state.  Since  measurements  are  made  at 
fixed  time  intervals,  the  Kalman  filter,  like  others  in  this 
report,  works  by  discrete  steps.  It  is  possible  then  to 
write  a state  transition  matrix  such  that  the  state  vector 
at  the  next  time  point  can  be  predicted  by  multiplying  the 
current  state  vector  by  a state  transition  matrix.  For  the 
example  system  this  is  nothing  more  than  a Taylor  series 
extrapolation . 


The  observed  quantities  can  be  predicted  from  the 
predicted  state,  i.e.,  the  observation  vector  is  a function 
of  the  state  vector.  The  actual  observations  are  measured 
and  the  differences  between  actual  and  observed  values  are 
multiplied  by  an  optimal  weighting  matrix  to  obtain  the 
corrections  which  will  update  or  improve  the  current  estimate 
of  the  state  vector. 

Obtaining  the  optimal  weighting  matrix  is  the  key 
to  the  Kalman  filter.  Briefly,  the  weights  are  derived  by 
consideration  of  the  covariance  matrices  of  the  observation 
and  state  vectors  so  that  that  state  vector  error  is  minimized 
in  a least-squares  sense. 

4.0  INDIVIDUAL  RANGE  FILTERS 

4.1  DET  1/SAMTEC  (ETR) 

4.1.1  Recursive  PV  Filter 

4. 1.1.1  Development  and  Usage 

A recursive  position  and  velocity  (PV)  filter  is 
the  first  filter  applied  to  the  input  raw  data  from  the 
radars  in  azimuth,  elevation  and  range.  This  filter,  called 
SMW,  is  used  to  pre-edit  the  raw  data  by  removing  wild 
points  from  the  input  sequence.  The  smoothed  position  and 
velocity  estimates  are  retained  merely  as  inputs  for  the 
next  call.  Data  received  by  the  computer  from  the  radars 
can  be  made  completely  erroneous  by  errors  in  transmission. 

If  such  data  were  allowed  into  the  velocity  filter,  the 
computer  velocity  would  be  bad  for  as  long  as  the  point 
stayed  in  the  filter  (normally  5 seconds).  This  in  turn 
would  cause  an  erroneous  impact  point  to  be  computed.  Since 
the  calculations  require  that  the  data  rate  be  maintained, 


the  edited  point  must  be  replaced  by  an  estimate.  For  this 
purpose  a light  linear  filter  seemed  most  desirable. 

4. 1.1. 2 Theory  and  Method  of  Operation 

A more  complete  discussion  of  the  theory  of  this 
form  of  filtering  can  be  found  in  Section  5,  Reference  13;  a 
synopsis  of  this  discussion  follows.  We  define  the  following 
parameters : 

{X(nT)}  = {Xfl}  = Raw  input  sampled  data  sequence. 

Sampling  is  at  equal  time  intervals  T. 

(W(nT)}  = {W  } = "Signal"  or  desired  component  of  the 

input  sequence  {XnK 

{y(nT)>  = {yn>  = "Noise"  or  undesired  component  of  the 

input  sequence  {Xn>. 

And  we  make  the  following  assumptions: 

a.  Let  X = W + y by.  definition 

n n n 

b.  The  expected  value  of  the  noise  sequence  yn  is 
zero,  and  its  autocorrelation  function  is  assumed  to  have 
the  form 

yAyj  = a2  if  i = j 
= 0 if  i f j 

c.  The  cross  correlation  between  the  signal  Wn  and 

the  noise  y is  assumed  to  be  zero, 
n 


k 
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It  is  desired  to  generate  a new  sequence  Pr  as  a function 

of  the  input  sequence  XR  by  means  of  a linear  time  invariant 

transformation  in  such  a way  that  the  new  sequence  is  a 

"best"  approximation  to  the  signal  WR  or  some  linear  function 

of  W . Hence.  P must  be  of  the  form 
n n 


P = 
n 


k=o 


gkXn-k 


( A-l ) 


which  is  a convolution  in  which  g^  is  the  desired  filter 
unit  impulse  response.  Alternatively  (A-l)  may  be  written 

as 


P(Z)  = G( Z )X( Z ) , 


( A-2  ) 


where  P(Z)  is  the  Z transform  of  PR , X(Z)  is  the  Z transform 

of  the  input  sequence  X and  G(Z)  is  the  Z transform  of  g, 

n k 

or  the  desired  filter  transfer  function. 


As  a form  for  a linear  transformation,  let 


g = Z a,  W . 
6nu  , k n-k 
k=o 


(A-3) 

be  an  arbitrarily  chosen  linear  function  of  the  signal 
which  it  is  desired  to  estimate  by  means  of  P , and  let  kQ 
be  the  ratio  of  steady-state  variance  in  PR  to  the  variance 
in  raw  input. 


The  variance  in  PR  is  given  by: 


r2  = “ 


= T - (Pn)2  = CT2 


(A-4) 


k=o 
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Therefore,  from  assumption  b, 


Let 


K 

o 


E 

k=o 


(A-5) 


D2  = ? (n  - E g (n-k))2  (A-6) 

n=o  k=o 

The  "best”  filter  is  defined  to  be  the  one  which  for  a 
given  Kq  will  make  D2  a minimum  or  for  a given  D2  will  make 
Kq  a minimum.  is  a measure  of  the  ability  of  the  filter 
to  reduce  the  random  noise  in  its  output.  D2  is  a measure 
of  the  ability  of  the  filter  to  respond  to  a signal  input, 
i.e.,  it  is  essentially  a measure  of  response. 


It  is  therefore  necessary  to  minimize: 


J = K + XD2 
o 


(A-7) 


where  A is  the  Lagrange  multiplier.  Substituting  the  values 
for  Kq  and  D2  in  (A-7)  yields 


I oo  oo 


J = 2 g2  + X 2 n-  2 gk(n-k) 


k=o 


n=o  n=o 


( A-8  ) 


The  optimum  position  filter  transfer  function,  G(Z),  is 
found  to  be 


£(Z)  - Z(qZ  + P-a)  . 

Z2-(2-a-/3)Z+(l-a) 


( A-9 ) 


The  velocity  filter  transfer  function,  G(Z)  , is  given  by 


G(Z)=  (L)  &SZ-W 


T Z2-(2-t*-0)Z+l-a) 


( A-10 ) 
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where 


« = 1 ZXZ2  (A-ll.l) 

S = (1-Zj.)  (1-Z2);  (A-11.2) 

and  Z^,  Z2  are  the  roots  of  the  quartic, 

(Z-l)4  + AT2Z2  = 0 , which  lie  inside  the  unit  circle 
in  the  Z plane.  (A-11.3) 

The  recursive  filter  set  of  equations  is  given  by 

(A-12.1) 

(A-12.2) 
(A-12.3) 
(A-12.4) 


Xn-XPn+a5n’ 

Xn  = Xn-l+^T5n' 
XPn+l=Xn+TXn 


To  edit,  6n  is  compared  in  absolute  value  with  a fixed 
edit  limit,  which  is  approximately  0.5  degree  in  angle  and 
400  feet  in  range.  If  it  exceeds  the  limit,  6n  is  replaced 
by  zero,  and  is  replaced  by  Xpn> 

4. 1.1. 3 Outstanding  Features 

Initialization 

This  filter  contains  a routine  for  self-initial- 
ization. Based  on  an  input  counter,  the  first  two  points 
are  differenced  for  the  initial  estimate  of  velocity  and 
the  second  point  is  used  for  the  initial  smoother)  position  estimate. 
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Consecutive  Edits 


If  five  consecutive  points  are  ediced,  the  filter 
resets  the  initialization  counter  to  call  for  self- 
initialization on  the  next  pass  through  the  filter. 

4. 1.1. 4 Limitations 


Since  this  filter  does  not  compute  acceleration, 
rapid  changes  in  velocity  will  cause  it  to  edit  and  reinitialize. 

4.1.2  Recursive  PVA  Filter 

4. 1.2.1  Development  and  Usage 

A recursive  position,  velocity,  and  acceleration 
(PVA)  filter  is  next  applied  to  the  edited  output  from  SMW. 

This  filter,  called  BSMW,  is  used  to  develop  estimates  of 
the  amount  of  noise  in  the  data.  These  noise  estimates, 
projected  to  impact,  are  later  used  as  a criterion  for 
selecting  the  radar  used  for  the  impact  calculations.  For 
the  purpose  of  obtaining  these  noise  estimates,  a function 
was  needed  which  could  reasonably  describe  the  motion  of  the 
vehicle  over  short  periods  of  time.  The  quadratic  was 
chosen  and  the  recursive  form  was  selected  for  savings  in 
time  and  core.  Light  filtering  is  used  in  order  to  minimize 
the  problems  associated  with  either  a staging  or  a malfunc- 
tioning missile. 


4. 1.2. 2 Theory  and  Method  of  Operation 


The  following  is  an  additional  synopsis  of  the 
discussion  found  in  Section  5,  Reference  13.  Using  the  same 
definitions  and  assumptions  given  for  the  PV  filter,  we  then 
let 


D2  = 


OO 

2 


OO 

2 


n=o\k= 


gkWn-k-8n5=2: 


OO 
2 

=o\k=o 


8kWn-k 


' 2 akWn-k 


k=o 


(A- 13) 
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The  function  to  be  minimized 


J,  becomes  in  this  case 


j-Sg’+X  S 2gkWn.l-fakW„-k| 
k=o  n=ok=o  k=o 


(A-14) 


The  optimum  PVA  filter  transfer  functions  are  then 


G(Z)  = Z|aZ2-  (2a#7)Z  + a *-r=j9|/(Z3-(3-a-0-7)Z2- 
(2a+0-r3)Z  - (1-a)] ; 

6(Z)  = Z(Z-l)  (0Z+27-0)/T[Z3-  (3-o-0-7)Z2  - 
(2a+/3-7-3)Z- (1-a)] 

G(Z)  = 27Z(Z-1)2/T2  [Z3  - (3-o-/3-7)Z2  - 
(2o+0-7-3)Z-  (l-a)]; 


(A-15.1) 

(A-15.2) 

(A-15.3) 


For  any  value  chosen  for  y,  the  roots  of  (Z-l) 6-AT4Z3=0 
which  lie  inside  the  unit  circle  in  the  Z plane,  namely  Z^ZgZ^, 
are  computed.  The  constants  a,g,y  may  be  computed  by  the 
relat ions : 


a - 1 - ZjZjZ^; 

0 = J4  [Z1Z2Z3  + -Z1-Z2-Z3-Z1Z2-Z3(Z1+Z2)] ; 

and 


( A-16 . 1 ) 
(A-16.2) 


7=1/2(l-Z1)(l-Z2)(l-Z3). 

(A-16. 3) 

Inversion  of  the  set  (A-15)  yields 

Iter  equations  in  the  time  domain: 

the  following  set  of 

«n  = Xn  - XPn; 

(A-17.1) 

Xn  = XPn+a6n 

(A-17.2) 

Xn  = Xn.1+TXn.1+0/T  6n; 

(A-17.3) 

Xn  = V1  +276n/T2; 

(A-17.4) 

XPn+1  =Xn+TXn+T2Xn/2 

(A-17.5) 
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This  filter  contains  an  edit  capability  of  the  same 
form  as  that  of  SMW,  but  since  a very  large  edit  limit  is 
specified,  this  capability  is  not  used. 

4. 1.2. 3 Outstanding  Features 

This  filter  also  contains  a routine  for  self- 
initialization. Based  on  an  input  counter,  the  first  three 
points  are  second-differenced  for  an  initial  acceleration 
estimate,  the  second  and  third  points  are  differenced  for 
velocity,  and  the  third  point  is  used  as  the  initial  smoothed 
position  estimate. 

4. 1.2. 4 Limitations 

Rapid  changes  in  acceleration  can  cause  this 
filter  to  either  load  or  lag  the  raw  data,  leading  to  an 
apparent  increase  in  the  noise  in  the  data. 

4.1.3  Finite  Span  Velocity  Filters 

4. 1.3.1  Development  and  Usage 

In  order  to  develop  a velocity  estimate  for  impact 
prediction,  the  data  from  the  selected  radar  is  passed 
through  one  of  three  filters:  GRSM46  (51  point,  half-second 
delayed,  for  10  pps  data);  GRSMPP  (51  point,  midpoint,  for 
10  pps  data);  or  G2SM00  (21  point,  midpoint,  for  20  pps 
data).  In  order  to  perform  this  task  in  a minimum  amount  of 
time,  it  was  desirable  to  take  advantage  of  a special  capa- 
bility of  the  machine  which  allows  the  binary  exponent  of  a 
floating  point  number  to  be  altered.  This  operation  effec- 
tively multiplies  or  divides  numbers  by  powers  of  two  at  a 
considerable  time  saving  over  floating  point  multiplication. 


20 


Therefore,  a method  which  consisted  of  a linear  combination 
of  data  points  utilizing  powers  of  two  was  sought  as  a 
substitute  for  a full  least-squares  solution.  The  finite 
span  method  was  chosen  over  the  recursive  in  order  to  get  a 
heavier  filter  without  the  overshoots  and  settling  problems 
inherent  in  the  recursive.  This  also  gives  a fixed  delay  in 
the  output,  whereas  the  recursive  is  generally  considered  to 
have  a variable  delay. 

4. 1.3.2  Theory  and  Method  of  Operation 

We  define  the  following  parameters: 

XR  = Input  data  sequence  sampled  at  equal  time 
intervals  T; 

Wr  = "Signal,"  or  desired  component  of  Xr ; 

Y = "Noise,"  or  undesired  component  of  X.. 
n n 

And  we  make  the  following  assumptions: 


Vk  ’ "n  + k™„  + *n'  <-NikiP>’ 

The  least-squares  fit  for  velocity  can  be  expressed  in  this 
form : 

k=P 

2 gkxn+k. 

k=-N 

These  three  filters  are  all  designed  to  be  used  with 
double-stored  tables  so  that  the  tables  do  not  need  shifting 
as  successive  points  are  added.  The  starting  location  of 
the  double-stored  table  and  the  index  to  the  oldest  point  in 
the  table  is  passed  to  the  routine. 
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4. 1.3. 3  Outstanding  Features 


These  three  filters  are  notable  fox  their  speed, 
both  in  required  overhead  and  actual  computation. 

4. 1.3.4  Limitations 

The  main  drawback  of  these  filters  is  the  storage 
required  for  the  tables  - 70  cells  for  each  coordinate  of 
the  double-stored  table. 

4.1.4  Variable  Edit  Filter 
4. 1.4.1  Development  and  Usage 

This  filter  was  developed  in  conjunction  with  the 
orbit  phase  real-time  program  at  ETR.  For  orbital  passes 
which  required  high  elevation  angles  of  the  radars  and 
consequently  rapid  changes  of  velocity  and  acceleration  in 
azimuth,  the  combination  of  SMW  and  BSMW  proved  to  have 
inherent  drawbacks.  It  seemed  desirable  then  to  transform 
the  data  first  to  a Cartesian  .coordinate  system  before 
attempting  to  edit  or  filter  it.  Once  the  data  were  trans- 
formed, the  next  problem  was  to  establish  a meaningful  edit 
limit.  For  short  ranges,  a given  angular  error  produced  a 
smaller  error  in  the  Cartesian  system  than  it  would  for  a 
longer  range.  There  are  also  effects  from  target  geometry 
which  cause  fluctuations  in  the  noise  level.  Consequently, 
a fixed  percentage  of  the  data  was  edited  by  allowing  the 
edit  level  to  vary.  It  became  apparent  that  this  edit  level 
was  then  a measure  of  the  noise  in  the  data  and  could  be 
used  as  a criterion  for  radar  selection.  Since  this  scheme 
replaced  the  SMW-BSMW  combination  and  several  routines  for 
error  averaging  and  transforming  to  the  Cartesian  system,  a 
considerable  time  savings  was  realized  as  a fortuitous  by- 
product . 


4. 1.4. 2 Theory  and  Method  of  Operation 


We  define  the  following  parameters: 

R = Ratio  of  editing  desired; 

E = Number  of  points  edited; 

P = Number  of  dat»  points  processed; 

L = Edit  limit. 

We  then  can  compute  the  following  quantities: 

C = E/P  = Ratio  of  editing  achieved  (B-l) 

M = C/R  = Edit-limit  multiplier  (B-2) 

The  smoothing  scheme  employed  is  the  same  as  that  of 
BSMW,  the  recursive  PVA  filter.  The  position  difference  is 
computed  as  in  equation  (A-17.1),  and  then  compared  in  the 
absolute  value  to  the  edit  lirtiit.  This  is  done  for  all 
three  position  coordinates. 

If  the  edit  limit  is  exceeded,  the  data  point  is  edited 
in  all  three  coordinates.  Once  the  point  is  edited,  the 
edit  limit  multiplier  is  checked  to  ensure  it  is  greater 
than  1.0  and  then  is  used  to  adjust  the  edit  limit. 

L = L • M (B-3) 

If  the  point  is  not  edited,  the  edit  limit  multiplier 
is  checked  to  ensure  it  is  less  than  1.0  and  then  is  used  to 
adjust  the  edit  limit  as  in  equation  (B-3). 
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4. 1.4. 3 Outstanding  Features 


In  order  to  keep  the  edit  limit  multiplier  M from 
becoming  unaffected  by  an  edit,  a maximum  is  set  upon  the 
value  to  which  P,  the  number  of  data  points  processed,  is 
allowed  to  go.  When  it  reaches  this  limit,  it  is  reset  to 
its  minimum  and  E,  the  number  of  points  edited,  is  adjusted 
so  that  their  ratio  C remains  essentially  the  same.  This 
process  maintains  a desired  level  of  "responsiveness”  in  the 
edit  level  itself.  If  several  points  in  a row  are  edited, 
the  edit  level  will  then  increase  in  an  almost  exponential 
manner  until  editing  ceases. 

4. 1.4. 4 Limitations 

This  filter  contains  no  internal  logic  for  reini- 
tialization. If  it  is  desired  to  reset  the  filter,  the  point 
counter  must  be  set  to  zero  externally. 

4.2  Space  and  Missile  Test  Center 

4.2.1  Exponential  filters- are  currently  used  in  the 
primary  flight  safety  computer  (IBM  7044)  and  the  backup 
flight  safety  computer  (CDC-3300).  The  same  basic  filter 
will  also  be  used  in  the  IBM  360-65  computer  when  it  replaces 
the  IBM  7044  as  the  primary  flight  safety  computer. 

4. 2. 1.1  Exponential  filters  were  incorporated  into  the 
Point  Arguello  Range  Safety  Impact  Predictor,  the  predecessor 
to  the  Vandenberg  Impact  Prediction  System  (VIPS),  in  June 
1965.  These  filters  were  then  carried  on  into  VIPS  when  the 
program  was  renamed  after  reprogramming  for  the  IBM  7044  in 
late  1966.  The  filtering  is  accomplished  in  two  subroutines 
in  VIPS.  They  are  EDIT  and  EDIT2 . EDIT  is  used  to  (1) 
compute  the  mean-average  deviation  between  the  predicted 
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value  and  the  measured  value,  (2)  edit  wild  data  points  from 
raw  data,  and  (3)  smooth  the  polar  coordinates.  EDIT2  is 
used  to  (1)  smooth  geocentric  position  and  (2)  differentiate 
this  position  to  obtain  geocentric  velocity. 

VIPS  is  the  primary  real-time  computer  program 
which  supports  all  major  launches  from  Vandenberg  AFB,  Cali- 
fornia. VIPS  operates  on  the  IBM  7044.  The  Backup  Impact 
Prediction  System  (BIPS)  is  the  secondary  real-time  computer 
program  which  operates  on  the  CDC-3300  computer.  BIPS  is 
equivalent  to  VIPS  except  for  automatic  abort,  number  of 
input  sensors  and  other  aspects.  The  filters  in  both  systems 
are  basically  identical.  Radar  measurements,  measurements 
from  the  General  Electric  Radio  Tracking  System  (GERTS)  and 
data  extracted  from  telemetered  inertial  guidance  systems  are 
used  as  inputs  into  VIPS.  These  data  are  checked  for  site 
identification,  parity  and  on-target  indicators  to  determine 
data  quality.  These  data  are  then  used  to  compute  impact 
prediction,  present  position  and  display  information  for  the 
Missile  Flight  Control  Officer  (MFCO). 

The  data  required  tb  compute  impact  prediction  for 
a missile  are  (1)  the  position  vector  and  (2)  the  .velocity 
vector  of  the  missile.  Three  independent  components  of 
position  and  three  independent  components  of  velocity  meet 
these  requirements.  These  data  are  obtained  by  direct  measure- 
ment from  the  GERTS  and  from  telemetered  inertial  guidance 
system  measurements.  GERTS  position  data  is  edited  to  remove 
wild  data  spikes  and  smoothed  to  reduce  the  noise  on  the 
measurements.  GERTS  velocity  measurements  are  used  in  raw 
form  without  either  smoothing  or  editing.  All  six  parameters 
from  telemetered  inertial  guidance  are  used  directly  and  are 
not  smoothed,  although  gross  editing  of  the  inertial  guidance 
prior  to  input  into  the  impact  prediction  routine  is  performed. 
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Most  radar  measurements  consist  of  the  three  compo- 
nents of  position  only  and  other  means  of  producing  velocity 
must  be  used.  Currently,  these  position  measurements  are 
edited,  smoothed  and  a cross-range  error  computed  for  use  in 
source  selection.  The  smooth  polar  position  components  are 
then  transformed  into  a geocentric  coordinate  system.  The 
geocentric  positions  are  then  smoothed  and  independently 
differentiated  to  produce  velocity  for  use  in  impact  prediction 
and  acceleration  for  use  in  prediction  for  updating  position 
and  velocity.  The  impact  prediction  is  then  corrected  for 
the  effects  of  drag  and  earth  oblateness.  The  final  predicted 
impact  point  is  then  displayed  for  the  MFCO. 

GERTS  data  (R,  A,  E,  R,  P,  Q)  are  input  and  the 
position  data  are  edited  and  smoothed  in  subroutine  EDIT. 

The  range  (Rq)  from  the  central  rate  station  to  the  object  is 
computed  and  edited. 

The  ranges  (R^  and  Rg)  to  the  other  two  rate  stations 
are  computed.  These  are  differenced  to  obtain  P and  Q. 


P = R - R1 
o 1 


Q - Ro  " R2 


These  quantities  are  then  edited  and  differentiated 
to  compute  R , P,  Q,  respectively  in  subroutine  EDIT.  If  the 
raw  input  range  rate  R is  not  good,  the  Rq  computed  in  EDIT 
is  substituted.  If  the  raw  lateral  rates  P and  Q are  good 
they  are  used  in  further  computations;  however,  if  the  raw  P 
and  Q are  not  good,  then  the  values  computed  in  EDIT  are 
substituted.  Geocentric  quasi-inertial  position  and  veloc- 
ities are  then  computed  to  be  used  in  impact  prediction 
computations . 
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4. 2. 1.2  The  inputs  for  EDIT  are: 

a.  Raw  data  argument  R azimuth  or  elevation  for 

radar . 

b.  Data  quality  Flag  - either  good  or  bad. 

c.  Initial  estimates  of  sigma  - (IMAD)  from  a 
priori  estimates. 

There  is  one  filter  for  each  of  the  raw  data  argu- 
ments, i.e.,  R,  A and  E for  each  radar  being  input. 

The  number  of  samples  of  filter  initialization 
(COUNT)  is  tested  against  an  optimum  number  (LIM1).  If  COUNT 
is  less,  the  edit  filter  undergoes  initialization.  During 
initialization,  COUNT  is  incremented  and  smoothing  constants 
a and  0 are  updated. 

COUNT  = COUNT  + 1 

a = cQ'jjflT  as  COUNT  approaches  LIM1 

6 = 1 - a 

Also,  if  COUNT  is  less  than  three,  no  mean-average 
deviation  updating  or  data  predicting  is  allowed.  If  COUNT 
equals  LIM1,  filter  initializing  is  bypassed  and  editing  is 
allowed.  If  COUNT  is  greater,  the  edit  filter  undergoes 
reinitialization.  During  reinitialization,  COUNT  is  reset  to 
LIM1  (value  - 20),  and  new  smoothing  operators  are  calculated: 


Si  ‘ Xn-I  ■ V„.!  Ip 


where 


COUNT  = LI Ml 

1_ 

“ COUNT 

6 = 1 - a 


If  COUNT  is  greater  than  three,  the  predicted  value  for  the 
current  sample  is  calculated: 


. = X„ 
pred  n 


, + V , + An-1 
■ 1 n- 1 


This  predicted  value  is  compared  with  the  raw  input 
value  and  a forecast  deviation  (DELXYZ)  is  calculated.  If 
this  forecasted  deviation  equals  zero,  editing  and  mean- 
average  deviation  updating  is  bypassed  and  data  smoothing  is 
performed.  Three  sigma  is  calculated  and  compared  with  the 
current  forecast  deviation.  If  the  forecasted  deviation  is 
greater,  the  raw  data  is  replaced  by  the  predicted  value.  If 
the  filter  is  being  initialized  (COUNT  is  less  than  LIM1), 
editing  is  bypassed  and  only  the  mean-average  deviation  is 
updated . 

a = + AXa  Mean-Average  Deviation 


where 


AX  = I X . - X I Forecasted  deviation  between 
1 pred  raw1 

predicted  and  raw  input  data. 

If  editing  is  performed,  the  forecasted  deviation 
is  compared  with  20  sigma.  If  the  deviation  is  greater  than 
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20  sigma,  the  mean-average  deviation  remains  unchanged  and 
the  edit  counter  (TALLY)  is  incremented.  After  the  raw  data 
is  thoroughly  evaluated,  data  smoothing  is  initiated.  New 
smoothing  operations  for  the  current  data  sample  are  calculated: 

si  - «*„ + »si-i 

Sn  = “S1  * 

Sn  ' “S?  ' 0SJ.1 


These  operations  are  linearly  combined  to  compute 
smoothed  position,  velocity  and  acceleration. 

X = 3(S*  - S2)  + S'*  smoothed  argument, 
n ' n n'  n 

Vn  = “ | (6- 5a)  S' -(10 -8a)  S2  + (4-3a)  ( smoothed  velocity 

2^2  of  argument 

An  = a2  (sj,  - 2S2  + smoothed  acceleration  of 

T 4P2~  " " argument 

The  smoothed  acceleration  will  then  be  used  in 

calculating  predicted  positiofi  (X  ,)  for  the  next  sample. 

preu 

EDIT2  - The  inputs  to  this  subroutine  consist  of: 

a.  Quasi-inertial  position  (XI,  ETA,  ZETA) 

b.  Filter  weight  factor  (a). 

The  outputs  consist  of: 

a.  Smoothed  quasi-inertial  position  (XIs,  ETAs , 

ZETAs ) 

b.  Quasi-inertial  velocity  (XIDOT,  ETADOT,  ZETADOT) 
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Quasi- inert ial  acceleration  components, 
a = Filter  weight 

P = 1 - a 


c . 


S 


* 

n 


s 


I 

n-1 


Argument 


= X 


n-1 


0Vl 


vn-1 


2<*Vn-1 

a 


4. 


Si  1 


^n-l 

a 


Data  is  smoothed  by  the  following: 


S2 

n 

s3 


Sn-l  ^ + S> 
Sn-1  ^ + ^ « 
Sn-1  ^ + S2  a 


First  operator 
Second  operator 
Third  operator 


Xn  = 3(sJ  - S2)  + Sj|  Smoothed  argument 

Vn  = |(8a-10)S2-(3a-4)S^-(5tt-6)sJj 


a 

7P  2 


Smoothed  velocity 
of  argument 


A = (S  - 2 S2  + s^i  rv2/  n 2 

11  n n id  ' p Smoothed  accelerator  of  argument 


4.2. 1.3  The  exponential  is  flexible  in  that  the  amount  of 
smoothing  can  be  changed  readily  to  minimize  noise  transmitted 
through  the  filter.  It  is  important,  however,  that  an  appro- 
priate compromise  be  made  between  how  much  noise  is  passed 
through  the  filter  and  how  fast  the  filter  is  able  to  respond 
to  the  data.  The  method  allows  for  (1)  efficient  data 
editing,  (2)  prediction  for  updating,  (3)  data  smoothing,  (4) 
data  differentiation,  and  (5)  an  evaluation  of  individual 
sensor  data  quality  for  use  in  source  selection. 


4.2. 1.4  The  limitations  of  this  type  of  filter  are  (1)  it 
is  fairly  difficult  to  initialize  and  (2)  it  is  not  responsive 
to  rapid  changes  in  acceleration. 
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4.3  KWAJALEIN  MISSILE  RANGE 


4.3.1  Kwajalein  Range  Safety  System  Filter 

The  Kwajalein  Range  Safety  System  (KRSS)  is 
designed  to  provide  a range  safety  capability  which  can  be 
applied  to  those  airborne  vehicle  test  programs  which 
require  range  safety  support  by  the  Kwajalein  Missile  Range 
( KMR ) . The  filter  developed  for  the  KRSS  is  presented 
herein.  The  characteristics  of  this  filter  are  described, 
and  the  method  adopted  for  selecting  the  optimum  operating 
point  is  given  in  detail. 

The  KRSS  filter  consists  of  an  expanding  memory 
polynomial  filter  of  degree  two  for  the  purposes  of  "state 
vector  initialization"  followed  by  a fading  memory  polynomial 
filter  of  degree  two  after  "state  vector  initialization." 
There  are  a number  of  reasons  for  choosing  this  filter 
design . 


a.  This  filter  combination  virtually  eliminates 
initialization  errors.  This  is  because  the  algorithm  is 
completely  self-starting.  After  precisely  3 cycles  the 
initial  values  will  have  been  completely  dropped. 


b.  There  is  an  excellent  theoretical  basis  for 
the  filter  equations  which  are  rigorously  derived  in  Section 
5.0,  Reference  (13).  Closed  form  expressions  for  the  variance 
reduction  factor  (VRF)  associated  with  position,  velocity, 
and  acceleration  are  available  for  both  the  expanding  and 
fading  memory  algorithms. 


c.  The  expanding  memory  polynomial  filter 
combined  with  the  fading  memory  polynomial  filter  can  be 
formulated  as  a set  of  extremely  compact  recursive  algorithms. 
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d.  The  Laguerre  discrete  orthogonal  polynomials 
form  the  basis  for  the  fading  memory  filter,  and  have  an 
exponential  weight  function  0 which  "fades"  out  in  a well 
behaved  manner  as  X increases.  This  parameter  0 is  useful 
in  controlling  the  balance  between  random  noise  suppression 
and  filter  dynamic  lag  and  in  characterizing  filter  behavior. 

e.  The  filter  design  facilitates  the  optimization 
of  the  filter  weights  during  real-time  operation  for  a 
variety  of  range  safety  applications. 

4.3.2  KRSS  Filter  Description 

The  KRSS  filter  equations  consist  of  an  expanding 
memory  polynomial  filter  of  degree  two  for  "state  vector 
initialization"  and  an  exponentially  fading  memory  polynomial 
filter  of  degree  two  after  "state  vector  initialization" 
has  been  completed.  The  equations  are  rigorously  derived 
in  Section  5.0,  Reference  (13)  and  are  structured  for 
recursive  use. 

4. 3. 2.1  Filter  Equations 

The  filter  equations  are  presented  in  the  following 
paragraphs : 

The  residual  is  given  by: 

ex(i)  = x(i)  - x(i/i-i);  x -*•  y,z  (1) 

The  current  state  vector  estimate  is  given  by  the 

following : 

x(i/i)  = x(i/i-l)+(2y/T2)ex(i)  (2) 

:*-( i/i ) = x(i/i-l )+(B/t  )ex(i)  x -►  y,z 

x(i/i)  = x(i/i-l)+(a)cx(i) 
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The  1-step  prediction  estimate  oi  the;  state 
is  given  by  the  following: 

x(  i+l/i ) = x(i/i ) 

x(i+l/i)  = x( i/i  )+( r )x( i/i ) x > y,z 

x( i+l/i)  = x(i/i)+(r)x(i/i)+(l?)x(i/i) 


(3) 


These  recursive  equations  give  the  filter  estimate 
of  the  state  vector  (X,  V,  A)  evaluated  for  the  current 
point  (Equation  2),  and  for  a 1-step  prediction  (Equation 
3),  are  valid  for  both  the  expanding  memory  and  the  fading 
memory  algorithms.  This  is  convenient,  since  the  only 
change  which  occurs  in  the  transition  from  expanding 
memory  to  fading  memory  filter  is  the  manner  in  which  the 
weighting  coefficients  (a,  6 and  y)  indicated  in  Equation 
(2)  are  calculated. 


For  expanding  memory  operation,  a,  3 and  y are  functionally 
related  to  the  number  of  cycles  of  data  processed,  N=0,l,2... 
For  fading  memory  operation,  a,  3 and  y are  functionally 
related  to  the  filter  control  parameter  8.  These  functional 
relationships  will  be  defined  below. 


The  following  definitions  apply  to  Equations  (1), 
(2)  and  (3): 


x(i)  = Measured  x position  at  time  t ^ (raw  inputs 
received  at  current  time). 

The  symbol  (~)  over  the  variable  denotes  that  the 
variable  is  an  estimate  derived  from  the  filtering  process. 
The  double  subscript  index  shown  in  Equations  (1),  (2)  and 
(3)  are  defined  as  follows:  the  first  index  refers  to  the 
instant  of  validity  of  the  variable  estimate.  The  second 


1 


33 


index  refers  to  most  recent  data  utilized  in  obtaining  the 
estimate.  For  example, 

x(i/i-l)  = Smooth  estimate  of  x position  for  time 
t^  obtained  from  data  taken  up  to  and 
including  that  at  time  t ^ ^ . 

x(i/i)  = Smooth  estimate  of  x position  for  time 
t^  obtained  from  data  taken  up  to  and 
including  that  at  time  ti- 

x(i+l/i)  = Predicted  estimate  of  x position  for 
time  ti+1  obtained  from  data  taken 
up  to  and  including  the  current  time 

t . . 

1 

The  velocity  and  acceleration  estimates  are  defined  in 
a similar  fashion,  utilizing  double  indexes  to  denote 
estimate  validity  time  and  data  validity  time. 

i = Time  interval  between  data  cycles. 

e (i)  = Residual  computed  during  current 

cycle . 

The  weight  coefficients  a,  6 and  y of  the  expanding 
memory  filter  are  functionally  related  to  N,  the  number  of 
cycles  of  data  processed,  by  the  following  equations: 


o(N)  = 

3(3N2+3N+2) 
(N+3)  (N+2)  (N+l) 

(position) 

(4) 

6(N)  = 

18(2N+1) 

(+3)  (N+2)  (N+l) 

(velocity) 

(5) 

y(N)  - 

30 

(N+3)  (N+2)  (N+l) 

(acceleration) 

(6) 
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With  N taking  on  successively  the  values  0,1, 2, 3,...  Nm-1, 

where  N is  the  value  of  the  counter  at  which  switching 
pi 

takes  place  from  the  expanding  memory  filter  to  the  fading 
memory  filter. 

The  weight  coefficients  a,  6 and  y of  the  fading 
memory  filter  are  functionally  related  to  the  fading  memory 
filter  control  parameter  0 by  the  following  equations: 


a(6)  = 1-03 

6(6)  = 3/2 (1-e)2 (1+9) 

Y(P)  = |(l-0)3 


0 < 0 < 1 


(7) 

(8) 
(9) 


4. 3. 2. 2 Variance  Reduction  Factor 

4 . 3 . 2 . 2 . 1 Expanding  Memory 

The  VRF  for  acceleration  for  an  expanding  memory 
filter  is  given  by: 


VRF  A = 


720 

t4(N+3)(N+2)(N+1 )N(N-1) 


(10 


This  factor  is  independent  of  where  the  evaluation  is  made 
(current,  1-step,  etc.).  This  is  because  of  the  choice  of 
a second  degree  polynomial  which  implicitly  assumes  that 
the  acceleration  being  estimated  is  constant.  Once  estimated 
for  a set  of  data  points,  its  estimate  and  the  associated 
VRF  is  independent  of  when  the  evaluation  is  made. 


The  VRF  for  velocity  is  dependent  upon  where  the 
evaluation  is  made.  For  a 1-step  prediction  this  is  given 
by: 

192N2+744N+684 

VRFV ( 1 ) = — o 

t2(N+3)(N+2)(N+1)N(N-1)  (11) 


The  "current  point"  VRF  for  velocity  has  a slightly  different 
form,  given  by: 


VRFV(O) 


192N2+24N-36 

t2(N+3)(N+2)(N+1)N(N-1) 


(12) 


Figure  1 shows  the  VRF  for  velocity  for  both  the 
current  point  and  for  a 1-step  prediction,  as  indicated  by 
Equations  (11)  and  (12). 


Finally,  the  VRF  for  position  is  even  more  depend- 
ent upon  where  the  evaluation  is  made  than  for  velocity. 

For  a 1-step  prediction,  the  position  VRF  is  given  by: 


VRFP(l) 


9N2+27N+24 

(N+l)N(Nrl) 


(13) 


and  the  current  state  VRF  for  position  is  given  by: 


VRFP(O) 


9N2+9N+6 
(N+l )N( N- 1 ) 


4. 3. 2. 2. 2 Fading  Memory 


(14) 


In  fading  memory  weight  Equations  (7),  (8)  and 
(9)  the  parameter  9 determines  the  weight  given  to  past 
data,  and  therefore  controls  the  tradeoff  between  dynamic 
response  and  noise  suppression  capability.  A large  value 
of  0 produces  good  noise  suppression,  but  poor  dynamic 
response,  whereas,  a small  0 gives  a good  dynamic  response, 
but  has  low  noise  suppression.  Figures  2,  3 and  4 show 
graphically  the  dependence  of  a,  8 and  y on  the  filter 
control  parameter  0. 
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FICUKE  1.  VELOCITY  VARIANCE  REDUCTION  FACTOR  VERSUS  "N"  FOR  EXPANDING  MEMORY  FILTER 


FIGURE  2.  ALl'HA  VERSUS  "T11ETA"  FOR  THE  FADING  MEMORY  FILTER 
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The  switching  point  N , between  the  expanding 
memory  filter  and  the  fading  memory  filter,  is  determined 
by  equating  their  velocity  variance  reduction  factors.  For 
the  second  degree  filter  selected,  this  point  is  given 
approximately  by: 

Nm  = j^. 463+(0-O.9)(4. 785-4. 463)/0. 099)J  /(1-0)  (15) 

Initialization  and  re-initialization  is  achieved 
with  the  use  of  the  expanding  memory  filter  weight  coeffici- 
ents and  a running  counter  N that  counts  successive  numbers 
of  "useable"  data  points. 


The  VRF  for  acceleration  for  the  fading  memory 
filter  is  given  by: 


VRF  A = 


6(1-8  )5 
t4( 1 + 0 )5 


(16) 


and,  as  with  the  expanding  memory,  it  is  independent  of 
where  the  evaluation  is  made. 


The  VRF  for  velocity  for  a fading  memory  filter 
1-step  prediction  is  given  by: 


VRFV ( 1 ) 


_1  (1-9 )3 
t2  (1+6)5 


( 49+500+1 302 

\ 2 


) 


(17) 


The  current  point  VRF  for  velocity  has  a slightly 
different  form,  given  by: 


VRFV ( 0 ) 


1_  (1-9)2  M9e2  + 5QQ  + i3 
t2  ( 1 + 9 )5  V 2 


(18) 


Figure  5 shows  graphically  the  functional  relation- 
ship given  by  Equations  (17)  and  (18). 
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Finally,  the  VRF  for  position  for  a 1-step  predic- 
tion is  given  by: 


VRFP(l)  = -1"9  , ( 19i-24e  + 1692+6e3+04 ) 

(1+0)°  ’ 


(19) 


and  the  VRF  for  position  evaluated  at  the  current  point  is 
given  by: 


VRFP(O)  = — ^ ( 1904  + 2403+1602+6e  + l ) 

(1-9)5 


(20) 


4. 3. 2. 3 Characteristics  Summary 

The  KRSS  filter  exhibits  the  following  character- 
istics : 


a.  The  filter  combination  is  recursive. 

b.  Starting  transients  will  be  essentially 
absent.  That  is,  transients  due  to  an  incorrect  a priori 
estimate  of  the  state  vector  will  be  eliminated. 

c.  The  algorithm  is  completely  self-starting. 
After  precisely  3 cycles  of  the  second  degree  algorithm, 
the  initial  values  will  have  been  completely  dropped. 

d.  At  each  cycling  of  the  expanding  memory 
filter,  its  error  covariance  matrix  corresponds  to  a least- 
squares  processing  of  all  previous  points.  At  each  cycling 
of  the  fading  memory  filter,  its  error  covariance  matrix 
corresponds  to  an  "exponentially  weighted"  least-squares 
fit  on  the  data  presented  to  it. 

e.  Dynamic  response  can  be  balanced  off  against 
random  noise  suppression  during  fading  memory  operation  by 
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the  appropriate  choice  of  the  filter  control  parameter  0. 
Setting  9 close  to  unity  is  best  from  a variance  reduction 
standpoint  and  worst  from  a dynamic  response  standpoint. 

The  reduction  of  0 from  unity  improves  dynamic  response  at 
the  expense  of  variance  reduction.  This  balance  can  be 
maintained  throughout  the  entire  tracking  interval  by 
appropriate  adjustments  to  0. 

f.  The  expanding  memory  filter  provides  a 
method  for  changing  in  a continuous  fashion  the  stress 
placed  on  the  data.  For  example,  a simple  process  for 
setting  the  stress  value  placed  on  reacquired  raw  data 
subsequent  to  a "data  dropout"  for  a period  of  time  can  be 
devised  so  that  the  stress  to  be  applied  is  determined  by 
the  length  of  the  dropout  interval  and  the  initialization 
state  at  the  time  dropout  occurred.  This  would  involve 
decrementing  the  data  cycle  counter  N by  appropriate  amounts 
during  the  dropout  interval  until  data  was  again  received. 

4.3.3  Filter  Operation 

4 . 3 . 3 . 1 General 

In  the  utilization  of  filtering  algorithms,  a 
compromise  must  be  made  between  the  conflicting  requirements 
of  good  noise  suppression  (heavy  filtering,  sluggish  system) 
and  of  responsiveness  to  changes  (light  filtering,  quick 
reacting  system).  This  compromise  applies  to  both  the 
expanding  memory  filter  and  the  fading  memory  filter.  For 
the  fading  memory  filter  with  a fixed  sample  interval  t the 
choice  of  the  filter  control  parameter  0 will  be  the  sole 
determinant  of  the  fading  memory  filter  response  to  noise 
and  systematic  changes  in  the  input  data  such  as  transients 
in  acceleration.  For  the  expanding  memory  filter  with  a 
fixed  sample  interval  x the  number  of  samples  N will  dictate 


44 


the  relative  weight  given  to  (1)  noise  suppression  and  (2) 
filter  responsiveness.  Initially,  for  small  values  of  N, 
the  expanding  memory  filter  will  be  extremely  responsive, 
with  correspondingly  poor  noise  suppression  characteristics. 

As  N increases,  filter  responsiveness  will  decrease  and 
noise  suppression  will  improve.  At  a specific  value  of  N, 
which  depends  on  the  fading  memory  filter  control  parameter  8 
(see  Equation  15),  the  transition  is  made  to  the  fading  memory 
filter.  Consequently,  8 is  the  single  parameter  which  can 
be  studied  to  specify  what  this  balance  must  be  for  a 
particular  application  for  the  expanding-fading  memory 
filter  combination. 

4. 3. 3.2  Filter  Performance  Indices 

Filtering  errors  are  caused  primarily  by  (1) 
filter  response  to  noise  on  the  input  data  and  (2)  the 
inherent  dynamic  lag  of  the  filter.  In  order  to  properly 
assess  these  two  independent  error  sources,  performance 
indices  have  been  adopted  which  depend  only  on  the  filter 
parameters  8 and  t;  where  8 is  the  fading  memory  control 
parameter  and  x is  the  input  data  interval. 


4. 3. 3. 2.1  Velocity  Variance  Reduction  Factor 

For  noise  suppression  the  generally  accepted 

performance  index  is  VRF.  For  range  safety  applications 

the  "velocity”  VRF  (rather  than  position  or  acceleration) 
% 

has  been  found  to  be  of  primary  importance  and  is  defined 
as  follows: 


VRFV(J)  = 


(J) 


_ Variance  in  Velocity  Estimate  at  Point 
Variance  in  Raw  Position  Input  Data 


and  has  the  units  of  (l/sec)‘ 


(21) 
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The  velocity  estimated  variance  is  dependent  upon 
the  point  being  evaluated.  If  the  velocity  variance  is  to 
be  evaluated  at  the  current  point  (J=0),  its  value  will  be 
lower  than  if  it  is  to  be  evaluated  for  a 1-step  prediction 
(J=l),  since  extrapolations  beyond  the  data  base  produce 
variance  increases. 

For  the  fading  memory  filter  this  factor  is  a 
function  of  J,  0 and  x denoted  functionally  by: 

2 

VRFV  ( J , 9,  T)  = ( J , 0 , i ) (22) 

°P 

to  emphasize  that  this  performance  index  is  an  explicit 
function  of  only  these  independent  variables.  However, 
subsequent  discussions,  the  simpler  form  (VRFV(J))  will 
used,  unless  functional  dependence  is  to  be  stressed. 

For  a 1-step  prediction  the  value  for  VRFV(l)  is 
given  by  Equation  (17).  For  a current  point  evaluation  for 
VRFV(O)  Equation  (18)  should  be  used.  This  current  point 
evaluation  of  the  velocity  variance  reduction  factor  has 
been  adopted  as  the  performance  index  which  characterizes 
the  noise  suppression  capability  of  the  KRSS  filter. 

Figure  5 shows  both  the  noise  suppression  perform- 
ance index  chosen  (VRFV(O)  = VRFV(0,  0,  x))  and  its  associated 
Ir-step  prediction  value  (VRFV(l)),  plotted  against  the 
filter  control  parameter  0.  This  functional  relationship 
will  be  used  subsequently  in  the  development  of  the  proper 
operating  value  for  0. 

4. 3. 3. 2. 2 Filter  Dynamic  Lag 

To  establish  an  index  of  performance  for  filter 


i n 
be 
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dynamic  lag,  a typical  systematic  change  in  the  input  data 
can  be  considered  to  be  a step  change  in  acceleration.  In 
the  absence  of  random  and  bias  errors  the  following  transient 
position  data  history  may  be  assumed  to  be  actual  raw  data 
for  each  component;  the  position  data  being  the  only  infor- 
mation acted  on  by  the  filter. 


x(  t ) 
x(  t ) 
x(t) 


0 

0 

0 


t>tA 


A (constant) 
A (t-tA) 

A ( t-t . )2/2 


Input  Change  In 

Acceleration 

Velocity 

Position 


Where  A is  the  acceleration  transient  step  magnitude  applied 
at  time  tA-  Emphasis  should  be  placed  here  on  the  fact 
that  A is  a step  or  change  in  acceleration  and  should  not 
be  confused  with  the  vehicle  total  acceleration. 

For  a particular  value  of  A,  a fixed  data  sample 

interval  t,  and  a choice  of  the  control  parameter  0,  the 

steady  state  (initialized)  fading  memory  filter  can  accept 

position  data  as  input  which  is  consistent  with  a step 

2 

input  in  acceleration  A,  i.e.,  x(t)=A(t-tA)  /2,  applied  at 
t=tA>  The  filter  response  to  just  this  step  transient  will 
exhibit  a velocity  error  (V^Aq)  which  will  eventually  damp 
out,  as  indicated  by  Figure  6.  The  time  delay  due  to  the 
filter  is  defined  as  the  time  interval  over  which 

a constant  acceleration  A acts  to  produce  the  maximum 
velocity  error  observed,  (a^lag^MAX‘ 


( Av  ) 

^ LAG ;MAX 


A (ATFILTER) 


However,  since  the  filter  is  linear,  the  filter  time  delay 
ATpiLTER  *s  independent  of  the  step  magnitude  A.  Functionally 
we  have: 


ATFILTER  ATFILTER(J’  9’  T) 


(23) 


-17 


VELOCITY 


1 


FIGURE  6.  ILLUSTRATION  DEPICTING  ELEMENTS  OF  FILTER  DYNAMIC  LAG 
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k J 


The  filter  dynamic  time  lag  (^FILTER 
upon  (1)  the  point  where  the  evaluation  is  made 
for  current  point,  J=1  for  a 1-step  prediction 
(2)  the  filter  control  parameter  9 and  (3)  the 
interval  x.  Therefore,  we  adopt  the  function 


) depends 
( i . e . , J=0 
evaluation ) 
sample 


atfilter(0)  ~ atfilter(0,  0’  T) 

as  the  performance  index  which  characterizes  the  dynamic 
lag  response  of  the  fading  memory  filter. 


Figure  7 shows  the  dynamic  time  lag  of  the  fading 
memory  filter  versus  0 There  is  a significant  difference 
between  the  lags  evaluated  at  the  current  point  (J=0)  or  at 
a 1-step  prediction,  with  the  lower  curve  (current  point) 
chosen  as  the  performance  index  of  KRSS  filter  dynamic  lag. 


4. 3. 3. 3 Velocity  Error  Minimization 


Defining  performance  indices  in  the  above  manner 
permits  the  calculation  of  a "combined"  velocity  error, 

AVFILTER’  as  follows: 

aVFILTER  = ^/a2(  aTf ILTER )^+  (N°p  VRFV 


where  ATpILTRR  and  VRFV  are 
performance  indices,  and  the 


the  previously  defined 
system  parameters 


f i Iter 


Op  = Standard  deviation  of  the  input  position  data 

N = Number  of  standard  deviations.  This  number 

indicates  the  weight  which  should  be  placed 
on  the  random  position  error  in  Equation  (24) 
When  the  step  acceleration  value  A used 
is  considered  to  be  a "maximum"  value,  a 
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0.1  SECOND 


FILTER  CONTROL  PARAMETER  THETA 


value  of  3 for  N is  typically  used.  If 
the  step  acceleration  value  A is  considered 
to  be  an  "average”  value,  a value  of  1 for 
N is  used. 

A = Magnitude  of  acceleration  step. 

For  a fixed  value  of  J (=0)  and  t (=0.1)  the 
value  of  the  filter  control  parameter  0 that  minimizes 
Equation  (24)  for  a specified  value  of  NOp  and  A is  assumed 
to  provide  the  best  balance  between  noise  suppression  and 
dynamic  response.  Equation  (24)  will  be  used  as  a basis 
for  selecting  the  filter  control  parameter  0 in  order  to 
optimize  the  filter  for  various  values  of  A and  Nap. 

Figures  8,  9,  10,  11,  12  and  13  illustrate  the 
utilization  of  Equation  (24)  to  calculate  the  velocity 
error  versus  8 for  a range  of  position  errors  (10,  50,  100, 
500  and  1000  meters)  for  six-step  acceleration  levels  of  1, 
10,  50,  100,  250  and  500  meters/second2,  respectively. 

Table  II  presents  a summary  of  the  minimum  values  of  velocity 
error  (AVMjN),  obtained  from  Figures  8-13  inclusive,  along 
with  the  associated  conditions  which  produced  each  result. 

4. 3. 3. 4 Control  Parameter  Selection 

The  optimum  value  of  the  filter  control  parameter 
0 is  that  value  which  minimizes  the  total  velocity  error 
given  by  Equation  (24).  Since  the  input  system  variables  A 
and  NOp  are  changing  as  functions  of  time  during  a mission, 
it  is  necessary  that  0 must  also  change  with  time  in  order 
that  the  optimum  value  be  continuously  selected.  In  fact, 
the  changes  which  occur  to  0 must  depend  upon  the  changes 
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50  METERS  100  METERS  500  METERS  1000  METERS 


POSITION  ERROR  50o  METERS  1000  METERS 


FILTER  CONTROL  PARAMETER  THETA 


POSITION  ERROR 
(N.  SIGMAP) 


FILTER  CONTROL  PARAMETER  THETA 


TABLE  II.  FILTER  PERFORMANCE 


d 


L 


A 

Nap 

A t 2 / N a p 

0 

avmin 

Refer 

to 

M/Sec2 

Meters 

- 

- 

M / Sec 

Figure 

1 

10 

0.001 

0.955 

2.2 

50 

0.0002 

0.975 

4.2 

100 

0.0001 

0.980 

5.6 

8 

500 

0.00002 

0.990 

10.7 

1000 

0.00001 

0.992 

15.7 

' 10 

10 

0.01 

0.885 

9.3 

50 

0.002 

0.940 

17.3 

100 

0.001 

0.955 

22.8 

9 

500 

0.0002 

0.975 

42.9 

1000 

0.0001 

0.980 

56 . 7 

50 

10 

0.05 

0.795 

25.4 

50 

0.01 

0.885 

46.4 

100 

0.005 

0.910 

60.6 

10 

500 

0.001 

0.955 

114.1 

1000 

0.0005 

0.965 

149.6 

100 

10 

0.1 

0.735 

39.4 

50 

0.02 

0.850 

71.3 

100 

0.01 

0.885 

92 . 8 

1 1 

500 

0.002 

0.940 

173.5 

1000 

0.001 

0.955 

228.2 

250 

10 

0.25 

0.635 

72.1 

50 

0.05 

0.795 

126 . 8 

100 

0.025 

0.840 

163.9 

12 

500 

0.005 

0.910 

303.2 

1000 

0.0025 

0.930 

397.5 

500 

10 

0.5 

0.550 

115.7 

50 

0. 1 

0.735 

197.0 

100 

0.05 

0.795 

253.7 

13 

1 

500 

0.01 

0.885 

464 .0 

1000 

0.005 

0.910 

606.4 
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which  occur  to  A and  Nop.  The  dependence  of  0 on  A and  Nop 
is  clarified  by  rewriting  Equation  (24)  in  the  following  form: 


' VFI LTER 


(25) 


The  nondimensional  parameter  [Ax2/Nop]  has  been 
found  to  characterize  the  effect  of  changes  in  A and  Nop  on 
the  minimization  of  AVpjpippR.  This  will  be  confirmed  from 
the  following  discussion.  The  quantities  [t2  VRFV ] and 
[ATPjjt£r/t1  can  be  considered  to  be  dimensionless  indices 
of  performance,  independent  of  x,  depending  only  on  the 
filter  control  parameter  0.  Equation  (18)  verifies  that 
[t2VRFV(0)]  is  indeed  independent  of  x.  The  dimensionless 
quantity  [ ATp TER/x 1 indicates  t^e  number  of  data  cycles 
of  dynamic  lag  resulting  from  a ' .ormalized"  step  accelera- 
tion A^2.  The  units  of  the  quantity  (Nop)  are  the  units  of 
velocity,  i.e.,  length/sec.  Consequently,  each  term  under 
the  radical  in  Equation  (25)  is  dimensionless,  and  independ- 
ent of  x.  A single  valued  functional  relationship  exists 
between  0 and  (Ax2/Nop).  That  is,  given  the  "input"  quantity 
(Ax2/Nop),  there  is  one  unique  value  of  0 to  minimize 
Equation  (25),  and  that  value  has  been  defined  as  the 
optimum  operating  point  for  the  fading  memory  filter.  This 
relationship  of  0 versus  (Ax2/Nop]  has  been  derived  for  the 
fading  memory  filter  and  is  shown  in  Figure  14.  The  rela- 
tionship given  in  this  figure  is  keyed  to  the  results  of 
Table  II. 

The  effect  of  0 on  noise  response  and  filter 
dynamic  lag  is  presented  in  Figures  5 and  7,  respectively. 
Using  Equation  (24)  and  the  relationships  represented  by 
Figures  5 and  7,  the  performance  curves  of  Figures  8 through 
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13  were  obtained.  For  each  case  there  is  a unique  value  of 
9 which  minimizes  the  velocity  error  of  the  filter.  These 
optimum  6 values  were  plotted  against  normalized  test 
conditions  in  Figure  14. 

Figure  14  shows  the  relationship  whereby  optimum 
filtering  can  be  maintained  continuously  during  a mission 
provided  detailed  information  on  the  time  (or  range)  depend- 
ence of  the  acceleration  transient  A and  input  position 

noise  error  Nc  are  known. 

P 

Figure  14  shows  that  optimum  filter  performance 
can  be  achieved  if  good  estimates  of  A and  Nop  can  be 
determined.  This  translates  filter  performance  parameters 
into  quantities  which  are  more  physically  related  to  the 
tracking  problem.  Generally,  op  can  be  estimated  as  a 
function  of  range.  The  value  of  A for  various  phases  of 
flight  can  also  be  determined  for  a particular  mission. 

The  choice  of  N,  the  weight  applied  to  the  position  error 
ap,  is  a constant,  subject  to  engineering  judgment.  When 
the  step  acceleration  value  A utilized  is  considered  to  be 
a "maximum"  value,  then  the  corresponding  value  of  3 for  N 
is  typically  chosen.  If  the  step  acceleration  value  A is 
an  "average"  value,  a value  of  1 for  N is  selected. 
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4.4  PACIFIC  MISSILE  TEST  CENTER 


4.4.1  Exponential  Filter 

4.4. 1.1  Development  and  Usage 

The  exponential  filter  as  used  at  the  Pacific 
Missile  Test  Center  (PMTC)  was  originally  developed  for  use 
on  ICBM/IRBM  class  vehicle  launches.  It  has  been  used 
subsequently  on  launches  of  shorter  range  surface-to-surface 
and  air-to-sur f ace  missiles.  The  filter  is  used  two  ways. 

It  is  used  on  raw  range  azimuth  and  elevation  data  from  each 
radar  to  provide  expected  values  of  the  next  observation  for 
use  in  editing.  The  mean  absolute  difference  between  the 
expected  and  observed  values  is  used  to  estimate  the  standard 
deviation  of  the  noise  as  a basis  for  editing  the  raw  data. 
The  noise  estimates  from  two  or  more  radars  are  then  compared 
to  determine  which  has  the  least  crossrange  Instantaneous 
Impact  Prediction  (IIP)  error.  This  radar  data  is  then 
selected  for  further  use. 

The  unsmoothed  but  edited  data  from  the  selected 
radar  is  converted  to  a Cartesian  coordinate  system.  A 
second  exponential  filter  operates  on  this  data  to  obtain 
smoothed  missile  position  and  velocity. 

4.4. 1.2  Theory  and  Method  of  Operation 

The  theory  of  exponential  smoothing  assumes  that 
the  impart  function  can  be  represented  by  a second  degree 

polynomial 

X.  = a + a.t  + a2  t2. 
t o 1 — o— 
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A future  value  can  be  predicted  by 


~ t n n 

L , „ = X.  + Z X + z X„. 
t+z  t t 2~  t 


If  a constant  (zeroth  degree)  model  is  assumed,  a recursive 
average  may  be  obtained  by 


Xt  = Xt-1  + a(Xt  " Xt-1)- 
where  a is  a smoothing  constant,  0<ct<l. 


Let  B = 1 - a 

and  sj  - it 

so  that  = B sj  j + aXt . 

S*  is  referred  to  as  the  first  order  smoothing  operator.  A 

C 2 
second  order  operator,  St , may  be  defined  by 

9 9 1 

St  B St-1  + “ St- 


Note  that  the  superscript  refers  to  the  order  of 
the  operator  and  is  not  an  exponent.  A third  order  operator, 

3 

St , is  similarly  defined. 

It  may  be  shown  that  linear  combinations  of  these 
three  operators  can  be  used  to  give  least-squares  estimates 
of  the  coefficients  of  the  polynomial. 


3SJ  - 3SJ  + s3t 


[ ( 6-5a ) sj  - 2(5-4a)S^  + (4-3a)SB] 


2iF  - 2Sl  + s?> 
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It  may  be  seen  then,  that 


X = A 

o 

X = Aj/At 

2 

X = A2/At 

4.4. 1.3  Outstanding  Features 

This  filter  provides  a simple  and  efficient  means 
to  fit  data  to  a second  degree  polynomial.  It  is  flexible 
in  that  the  value  of  a may  be  selected  to  give  the  desired 
compromise  between  response  and  smoothness. 

4.4. 1.4  Limitations 

The  filter  is  general  and  assumes  nothing  other 
than  the  fit  to  a second  degree  polynomial.  Additional 
smoothness  or  response  can  be  obtained  by  more  specialized 
f ilters . 

4.4.2  Recursive  Exponentially  Weighted  Least-Squares 

Filter 

4 . 4 . 2 . 1 Usage 

The  filter  is  exponential  and  acts  as  a one-step 

\ 

predictor.  It  is  an  expanding  memory  filter  until  the 
number  of  cycles  reach  the  filter  size  number  being  used  at 
that  time.  The  filter  switches  from  an  expanding  memory  to  a 
fading  memory  when  the  number  of  cycles  is  equal  to  or 
greater  than  the  filter  size.  The  telemetered  vehicle's 
angular  rate  information  is  used  to  control  the  filter  size 
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used  during  powered  flight.  For  example,  if  a bad  angular 
rate  signal  is  received,  then  the  short  filter  size  is  used, 
making  the  filter  output  noisier,  but  more  responsive  than 
when  the  long  filter  size  is  used. 


4. 4. 2. 2 Theory  and  Method  of  Operation 

The  theory  of  exponential  smoothing  assumes  that 
the  impact  function  can  be  represented  by  a second  degree 
polynomial . 


= ccq  + a^t  + dgt  / 2 


A future  value  can  be  predicted  by 

pt+z  - pt  + zvt  + z2V2 

The  equations  used  are  obtained  by  filtering  a 
second  degree  polynomial  in  time  to  the  position,  velocity, 
and  acceleration  data  using  a recursive  exponentially  weighted 
least-squares  procedure.  The  equations  are: 


APi  " °i  - Ai-1 


Ai  = [Ai-l  + VPi]K 


Qi  = F.C  + A. 


Vi  " Vi-1  + 2Qi  + 6iAPi 


Pi  = Pi-1  + Vi  - Qi  + aiAPi 


where 


= Observed  edited  Cartesian  coordinate  at  time 


obtained  from  data  taken  up  to  and  including 
that  at  time  t ^ ^ . 


N 


and 


= Smooth  estimate  of  velocity  for  time  t^ 

obtained  from  data  taken  up  to  and  including 
that  at  time  t^_^;  the  first  time  derivative 
of  position. 

= Smooth  estimate  of  acceleration  reflecting 

the  thrust  of  the  vehicle  for  time  t.  obtained 

1 

from  data  taken  up  to  and  including  that  at 
time  t^_^;  the  second  time  derivative  of 
position . 

F^  = Smooth  estimate  of  acceleration  of  the  vehicle 
due  to  earth's  gravity  for  time  t.  obtained 
from  data  taken  up  to  and  including  that  at 
time  t^ 

C = Conversion  factor  to  make  the  acceleration- 
due-to-gravity  term  dimensionless 

' 

K = Forcing  term  used  to  minimize  filter  overshoot 
during  the  main  engine  cutoff  (MECO)  phase  of 
the  flight;  this  is  a parameter  entry  whose 
value  is  approximately  unity. 

cij  = 3(3J2  + 3J  + 2 ) / ( J + 3 ) ( J + 2 ) ( J + 1) 

6j  = 18(2J  + 1)/(J  + 3 ) ( J + 2 ) ( J + 1) 

Yj  = 30/ ( J + 3 ) ( J + 2 ) ( J + 1) 
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where 


a,  8 and  y are  the  weight  coefficients  where 

J = 0,1,2 N.  N is  the  predetermined  filter  size  entered 

as  a parameter  before  the  operation.  N may  have  two  different 
values  during  powered  flight;  one  value  signifying  a long- 
filter  length  for  good-quality  telemetry  data,  a second 
value  signifying  a short-filter  length  for  bad-quality 
telemetry  data  and  a third  value  (usually  extremely  large) 
for  use  during  ballistic  flight. 

The  following  coefficients  are  used  during  ballistic 

flight : 


ctj  = 2(2J  + 1)/(J  + 2 ) ( J + 1) 

8j  = 6/ ( J + 2 ) ( J + 1) 

The  coefficient  y is  not  used  since  the  thrust  is  zero. 

4. 4. 2. 3 Outstanding  Features 

The  position  data  are  rotated  into  an  orientation 
such  that  the  filter  efficiency  is  maximized  by  minimizing 
the  random  and  systematic  radar  tracking  errors.  This 
orientation  is  in  the  theoretical  direction  of  the  vehicle 
at  MECO  from  the  radar. 

The  filter  takes  advantage  of  the  telemetry  infor- 
mation during  flight  and  utilizes  the  time  of  MECO,  which  is 
determinable.  The  quality  of  the  signal  is  determined  by  the 
values  of  the  pitch  and  yaw  rates.  The  filter  may  also 
minimize  overshoot  error  by  eliminating  residual  thrust 
after  burnout  is  sensed. 


4. 4. 2. 4 


Limitations 


This  filter  is  designed  to  be  used  for  vehicles 
going  through  a fairly  continuous  powered  flight  which, 
after  becoming  ballistic,  remains  so  throughout  the  remainder 
of  flight.  With  modifications,  however,  the  capability  of  a 
reignition  stage  with  powered  flight  followed  by  ballistic 
flight  can  be  incorporated.  Much  analysis  must  be  done  on 
the  vehicle  in  flight  so  that  the  many  parameters  can  be 
computed.  This  filter  is  not  as  generalized  as  the  one 
previously  described,  but  it  may  present  a more  true  estima- 
tion of  the  vehicle  behavior. 

The  filter  does  need  the  telemetry  data  in  conjunc- 
tion with  the  radar  data  to  be  smoothed.  Certain  parameters 
are  necessary,  such  as  the  orientation  angles,  for  error 
minimization  and  approximate  time  of  burnout  as  well  as 
vernier  engine  thrust. 

4.5  WALLOPS  ISLAND 

4.5.1  Free-Flight  Filters 

These  filters  were  designed  and  reported  by  the 

Special  Information  Products  Department  (SIPD),  Defense 

Electronics  Division,  General  Electric  Co.,  Syracuse,  New 

York.  They  are  an  integral  part  of  the  Instantaneous  Impact 

Prediction  System  supplied  by  General  Electric  in  1965  under 
% 

contract  NAS  6-852.  General  Electric  designed  the  filters 
around  Wallops  radar  data.  Since  its  development  minor 
changes  have  been  made,  but  the  basic  filter  has  been  pre- 
served. These  filters  utilize  the  method  developed  by 
Rudolf  Kalman.  These  filters  are,  however,  special  adaptations 
of  the  Kalman  method  in  order  to  meet  real-time  requirements. 
For  a detailed  explanation  of  this  technique  refer  to  refer- 
ences 1 through  6 of  Section  5.0. 


1 
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4. 5. 1.1  This  filter  was  designed  and  is  used  for  extended 
periods  of  free  flight  or  vehicle  coast. 


[ 


4.5. 1.2 


We  define  the  state  vector  of  the  system  to  be: 


X 


x 

y 

z 

X 

y 

z 


where  each  element  is  in  the  inertial  coordinate  frame. 

In  order  to  predict  the  state  from  one  data  sample  to  the 
next,  we  use  the  matrix  equation 

X( n/n-1 ) = !(n,n-l)  X (n-l/n-1)  + T W (n-1)  (2) 


where  W is  the  gravity  model  which  is  covered  in  subparagraph 
4.5.3.  The  notation  A(n/m)  will  be  used  to  designate  an 
estimation  of  A at  time  n based  on  observations  through  time 
m. 


For  our  system  we  know  the  following  scalar  equations: 

X(n/n-l ) = x(n-l/n-l ) + ti(n  = l/n-l)  + ll  x (n-1)  } 
x(n/n-l)  = x(n-l/n-l )+tx(n-l ) } 

where  t is  the  time  difference  between  data  samples  and  x 
denotes  a time  derivative.  Similar  equations  can  be  written 
for  y and  z.  These  scalar  equations  can  be  put  into  the 
matrix  form  of  Equation  (2)  to  yield  the  <f>  and  T matrices. 
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The  state  vector  is  related  to  the  measured  data  by  the 
equation 

Z(t)  = HX( t ) + V ( t ) (5) 

where  Z(t)  is  the  measurement,  H is  the  transformation 
matrix  and  V(t)  is  the  measurement  error  vector.  Since  we 
are  measuring  range  (R),  azimuth  (A)  and  elevation  (E),  this 
equation  is: 


Since  our  state  vector  is  in  the  inertial  frame,  we  must 
transform  from  there  to  fixed  earth,  to  radar  Cartesian  and 
finally  to  radar  coordinates. 
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r 


Let  4^  - Radar  longitude  from  Greenwich 


tp  = Program  time 
r = Earth's  rotation  rate 


* = Og  + rtp 


Then  we  have 


w 


cos  <p  sin  0 0 

-sin</>  cos  <f>  0 

0 0 1 


X| 

X 

1 

e 

>1 

-Aj 

0 

_zl_ 

(7) 


(n/n-1) 


(n/n-1) 

This  transformation,  including  a derivation  of  A,  is  discussed 
in  detail  in  subparagraph  4.5.4. 


Since 


R( n/n-1 ) = [u2(n/n-l ) + v2(n/n-l)  + w2(n/n-l)]“ 

- ‘“-1  [%mi\ 

^w( n/n-1 ) j 


(8) 


E(n/n-l)  = tan 


where  p = [u2(n/n-l)  + v2(n/n-l)]^ 


Since 


Ru 


dR 

nr 

dR 

17 

dR 

dw 

dA 

dA 

dA 

du 

dv 

dw 

dE 

dE 

dE 

du 

dv 

dw 

(9) 
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This  becomes 


R » 


U"  \ u p 

Rjp  R-p  ~r- 


Ii  - l<j  A 


CCS0  sin0  0 


-S1H0  co>0  0 

o 0 1 


0 0 0 
where  0=0  0 0 


Ip  0 o 


"?  ""  ff°rmUlate  the  fUter  e<,Uatl°"s-  ">•  "Iter  consists 
Of  the  following  three  equations: 


-(n/n-D  - X(n/n-l ) + K(n)  AZ(n) 

K(n)  - P(n/n-l)HT  |(g(n/n-l )_HT  * s_.  |-> 
P(n/n~  l ) = »,r„T  |ig-HT  + ? g ^ 
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where 


= Covariance  matrix  of  the  state  vector 

= P( n- 1 / n-2  ) 

= Difference  between  measured  (R,A,E)  and 
predicted  (R,A,E) 

= Measurement  covariance  matrix 

= Weighting  matrix 

= Covariance  of  W (the  gravity  model) 

which  we  assume  to  be  zero  (i.e.,  perfect 
gravity  model). 

Equation  (15)  reduces  to 

P(n/ii-l ) - <t>P  <1>!  <I>P  H 1 |HP  HT  + S J 1 HP4>T 

or 

R(n/n-1 ) = <t>p  >i> 1 ‘«>KvnjHP 

Now  define 
Mu  IHP  Ht  f s2|  1 

Then 

K(n)  = P hT  Mu 


(16) 


(17) 


(18) 


P 

P~ 

AZ(  n ) 

S 

— z 

K(n) 

Q 
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Now,  referring  to  Equation  (13),  we  have  the  best  estimate 
of  the  state  given  all  data  up  to  "now": 


Z'  n/n)  = _X(n/n-l ) + &fn) 


R ( n ) - R(n/n-l ) 

A(n)  - A(n/n-l ) 
E'n)  - E(n/n-l) 


(19) 


Where  X(n/n-l)  is  from  Equation  (4). 

R(n),  A(n),  E(n)  are  corrected  rate  measurements  now. 

R(n/n-l),  A(n/n-l),  E(n/n-l)  are  from  Equation  (8j  and  K(n) 
is  from  Equation  (18). 


As  it  is  used  by  NASA  at  Wallops  Island,  Virginia, 
the  free-flight  filter  depends  on  the  powered-flight  filter 
(subparagraph  4.5.2)  for  initialization  of  the  state  vector. 
An  initial  estimate  of  the  state  vector  is  computed  from  the 
output  of  the  last  cycle  of  the  powered-flight  filter.  The 
covariance  matrix  is  initialized  as  follows: 


106  0 
0 106 
0 0 
0 0 
0 0 
0 0 


0 0 
0 0 
106  0 
0 104 
0 0 
0 0 


0 0 
0 0 
0 0 
0 0 
104  0 
0 104 


4.5. 1.3  Outstanding  Features 

a.  Both  the  free-flight  and  the  powered-flight 
filters  are  end-point  filters,  which  is  essential  for  real- 
time applications. 
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b.  Both  filters  are  recursive;  thus,  only  one 
data  sample  must  be  stored  in  the  computer  at  any  one  time. 

4. 5. 1.4  Limitations 

Both  the  free-flight  and  the  powered-flight  filters 
are  sensitive  to  the  data  editing  scheme. 

4.5.2  Powered-Flight  Filter 

4. 5.2.1  This  is  the  primary  range  safety  filter  since  it 
is  used  during  the  powered  portion  of  flight  and  for  short 
periods  of  coast. 

As  burnout  is  approached  during  a powered  stage 
the  uncertainty  of  the  acceleration  terms  in  the  powered 
flight  model  is  increased.  This  is  done  by  setting  the 
variance  of  R , A and  E to  a relatively  large  value.  This 
action  enables  the  filter  to  "turn  the  corner"  in  its  calcula- 
tion of  velocity.  The  uncertainty  is  kept  large  until  after 
the  next  stage  has  fired.  Then  the  filter  covariance  matrix 
is  again  allowed  to  converge  until  the  burnout  of  this  stage 
is  approached.  The  reason  for  the  covariance  being  kept 
large  during  free  flight  is  the  immediate  response  of  the 
filter  to  accelerations  during  this  period.  This  may  be  due 
to  programmed  or  premature  ignition. 

NOTE;  For  the  5-second  interval  from  T-2.5  sec.  to  T+2 . 5 
% 

sec.,  where  T is  a nominal  staging  time, 

PR  = 1000.  x PR 
K33  .33 

P = 100.  x P 
33  33 

P„  = 100.  x Pr 
E33  E33 
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where  PR , PA  and  PE  are  the  covariance  matrices  for  range, 
azimuth  and  elevation,  respectively.  ^or  example, 


°*R  °Rd  °Rb  °Ru  °l 


R ReR 


4. 5. 2. 2 This  recursive  filter  uses  the  third  degree  polynomial, 


R = R0  t Rt  + ft—  * Rll 


where  t is  time  and  the  "dot"  stands  for  time  differentiation 


and  the  second  degree  polynomials 


" t ^ 

A = A + At  + A75— 
o 2 


E = E + Et  + E^ 
o 2 


for  azimuth  ar.d  elevation  models,  respectively.  To  prevent  a 
bias  in  the  outputs  of  the  filters,  data  is  weighted  out  in 
exponential  fashion  and  an  error  term  is  introduced  into  the 
model.  The  state  vectors  then  become 


R 

R 

R ,A  = 
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where  cR,  r:ft  and  are  the  error  terms  in  range,  azimuth 
and  elevation,  respectively. 

The  following  will  apply  only  to  range,  but  similar  arguments 
apply  to  azimuth  and  elevation.  We  would  like  to  write 
Equation  (20)  in  the  form  of  the  matrix  equation 

R(n/n-l)  = £(n,  n-l)R(n-l/n-l)  + W(n)  (24) 


where  the  notation  A(n/m)  refers  to  an  estimate  of  some 
quantity  A at  a time  n based  on  observations  up  to  and 
including  time  m.  The  W(n)  is  the  system  forcing  function 
which  we  assume  to  be  zero  with  a covariance  of  zero.  We 
can  write: 


R(n/n-f) 


0 


R(n-l/n-l)  W 


where  aR  is  a constant  (presently  equal  to  0.5).  The  state 
vector  is  related  to  the  measured  data  by  the  equation 


Z(t)  = HX( t ) + V(t) 


(201 


where  Z(t)  is  the  measurement,  H is  the  transformation 
matrix  and  V(t)  is  the  measurement  error  vector.  Since  we 
are  measuring  range,  we  have 
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r 


r 

R = (1000  01  j R 


so  a = (i  o o o o|. 


The  filter  consists  of  the  following  equations: 


R(n/n)  = R ( n / n - 1 ) + K(n)  |R(n)  - R(n/n-l)] 
K(n)  = P( n/n-1 ) HT  | H P ( n / n - 1 ) HT  + SR]_ 
P ( n / n - 1 ) = ^{P"  - P"  HT  I HP  HT  + SrJ_1.HP"}0 


where 


Covariance  matrix  for  range  state 


vector 


P(n-l/n-2) 


R(n) 


K(n) 


Present  range  measurement 


Variance  of  measured  range 


Weighting  matrix 


R(n/N-1)  is  computed  from  Equation  (25) 


The  filter  is  initialized  as  follows: 


£(1/0) 


10S 

0 

0 

0 

0 


R 


MO) 


0 

10,000 

0 

0 

0 


0 

0 

1000 

0 

0 


0 

0 

0 

1000 

0 


0 

0 

0 

0 

1000 


(31) 


where  Rq  is  the  first  range  measurment . 


If  the  filter  is  initialized  at  some  point  other 
than  during  the  first  three  cycles,  the  covariance  matrix  is 
initialized  as  above,  but  ft  and  R in  the  state  vector  are 
computed  from  first  and  second  differences  of  range  measure- 
ments . 


The  state  vector  R was  used  in  the  above  discussion, 
but  the  formulation  remains  the  same  for  azimuth  and  elevation 
with  the  following  exception  for  the  initialization  of  the 
covariance  matrix: 


£(i/0)  = 


10SA  ooo 

0 10  0 

o 0 .1  o 

o o o .1 


where  SA  is  replaced  by  S£  for  the  elevation  covariance 
matrix. 
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4. 5. 2. 3 Outstanding  Features 
Refer  to  subparagraph  4.5. 1.3. 

4. 5.2.4  Limitations 
Refer  to  subparagraph  4.5. 1.4. 

4.5.3  Acceleration  Due  to  Gravity 

Assuming  that  the  earth  is  a pear-shaped  oblate 
spheroid,  the  gravity  potential  (Fisher  Potential  Model)  can 
be  expressed  as  a function  of  geocentric  co-latitude,  the 
complement  of  geocentric  latitude  and  the  radial  distance 
from  the  center  of  the  earth. 


The  gravity  potential  formula  is: 
v(r-  *'c)  = TT  ' + | cos2  <t>'c  - ^ j 


DR- 


35  r- 


35  cos4  <t>' 


30  cos2  <J>|.  + 3 


H 


1 

■ 

cos^  4>' 

3 

COS  <I> 

l c 

5 c 

Earth  gravitation  constant 
1st  oblateness  correction 
2d  oblateness  correction 
Pear  shapness  correction 
Geocentric  radius 
Equatorial  radius  of  the  earth 

Geocentric  Colatitude 


Axis  in  equatorial 
plane  of  the  earth 
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The  gravity  acceleration  components  can  be  obtained  by 
taking  the  negative  gradient  of  the  potential  V(r,z),  that 
is, 

r = f (X,  Y,  Z) 
then  V( r , Z ) = U|f(x,y,z)z] 


and 


then 


and 


dV 

IT 


+ 


8V 

8z 


Of  _ 
9\ 


GM 

Rc 


= (>M 


2JR^Z 


_ v 

3f 

r 

• Tz  ~ 

+ j r3 

[-f- 

6Z2 
r 8 

+ 7r* 

+ DR5  | 

[ 4z3 

£ 

r 


4.5.4 

4.5.4. 1 
Position 


Transformation  Routines 

Transformation  from  Radar  Measurement  to  Rectangular 
Coordinates 


From  Figure  15,  the  transformation  from  radar  measurements 
R,  A and  E to  rectangular  position  coordinates  u,  v and  w is 


u - R Cos  E Sin  A East  positive 

v - R Cos  E Cos  A North  Positive 

w - R Sin  E'  Up  positive 


Tangent  plane  radar  oriented. 


(32) 
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Radar 


Figure  15.  Radar  Coordinate  System. 

A small  change  in  any  one  of  the  variables,  u for  example, 
could  be  expressed  in  terms  of  a Taylor  series  expansion 
about  the  mean  value  u for  small  deviations  as 


“Is- 


|_  + AE  i*iL 

A 9E 


when  higher  order  terms  are  neglected.  This  expansion  for 
all  three  variables  u,  v and  w may  be  written  in  matrix  form 


Designating  the  matrices  in  Equation  (34)  as  AU  = UUAR,  we 

K 

can  square  the  equation  as  follows: 


Au  AuT  = UR  AR(URAR)T  = UR  AR  ArT  UrT 

The  ensemble'  average  of  Equation  (35)  yields  a covariance 
matrix  of  the  u,  v and  w variables  in  terms  of  the  R,  A and 
E covariance  matrix  and  the  matrix  of  partials.  The  expan- 
sion of  Equation  (35)  is  carried  out  below. 

Squaring 


A u AuT  = UrAr(UrAr)t  = UrArArtUrt 


Au  Av  Au  Aw 


Av  Au  Av^ 


Aw  Au  AwAv  Aw  ^ 


Taking  the  ensemble  average 


0uv  ouw 


av^  Ovw 


owv  aw 


S1  = URSRUR1 


ARAE  ARAA 


Av  Aw  = UR  AEAR  AE2  AEAA  Ur 


AAAR  AAAE  AA' 


0 UR 1 Ali  cross 

correlations 
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Normally  Equation  (36)  is  written  in  the  general  terms 


S = M S M 
y xxx 


where 


Mx  is  the  matrix  of  partials. 


This  form  will  be  used  throughout  the  remainder  of  the 
report . 


The  elements  of  the  matrix  of  partials  M from  Equation  (32) 
are : 


uD  = cos  E sin  A 
K 


u^  = R cos  E cos  A 


uT,  = -R  sin  E sin  A 
E 


vR  = cos  E cos  A 


= R cos  E sin  A 


v„  = -R  sin  E cos  A 
E 


wR  = sin  E 


wA  = 0 


wT,  = R cos  E 
E 
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\ 


(37) 


f 38  i 


The  errors  in  the  u,  v and  w coordinates  are 


o 2 
u 

= 

2 2 9 9 9 9 

UR  aR~  + uE~aE  4 UA~  aA~ 

0 * 
uv 

= 

VR  UR  aR~  + 'eUE°E2  + vAua  °A2 

0 

uw 

= 

WRUR°R2  + WEUE°E2 

9 

0 - 
V 

= 

9 9 9 9 9 9 

VR~  aR~  + VE  °E  + vA“aA“ 

(40) 

o 

wv 

= 

WR  VR  aR“  + VE  WE  °E“ 

o 2 
w 

= 

9 9 9 9 

Wr  or“  + wE-  aE“ 

These  are  the  uncertainties  in  the  calculated  u,  v and  w set 

from  the  measured  R,  A and  E coordinate  set.  Because  the 

covariance  matrix  is  symmetrical,  o , a and  a are  not 

vu ’ wu  wv 

calculated. 


4. 5. 4. 2 Transformation  of  Position  from  Radar  to  Earth 
Centered  Fixed  Coordinates 


The  transformation  of  coordinates  from  u,  v and  w to  x,  y 
and  z may  be  expressed  as  three  rotations.  The  angles  of 
rotation  are  shown  in  Figures  16  and  17.  Performing  the 
transformations,  we  readily  obtain  the  following  results: 


1.  First  rotation  (Figure  17a) 


cos  a 
-sin  a 

0 


2. 


Second  rotation 


cos  6 w 0 

0 1 

sin  6w  0 


U* 

V' 

W’ 


(Figure  17b) 


-sin  6 w 

u 

0 

V 

cos  6 w 

w 

m 

» — 

(41  ) 


(42) 
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3. 


Third  rotation  (Figure  17c) 


0 -sin 

1 0 

0 COS0J 


cos  0^ 

U' 

0 

V' 

sin  0j 

w' 

(43) 


where  = 9*  + 6 . By  substituting  Equations  (41)  and  (42) 
into  Equation  (43),  the  following  results  are  obtained: 


X 


1“  — 1 

r* 

— 1 

r— 

—1 

r 1 

X 

0 

-sin  0 j 

COS  0] 

cos  6 w 

0 

-sin  6 w 

cos a 

sin  a 

0 

u 

y 

= 

i 

0 

0 

0 

1 

0 

-sin  a 

cos  a 

0 

V 

Z 

0 

COS  0 J 

sin  0 j 

sin  6 w 

0 

cos  6 w 

0 

0 

1 

w 

-J 

L 

J 

L J 

(44) 


or  by  multiplying  the  three  matrices  above;  the  result  is 
shown  in  Equation  (45). 


X 

rd 

r1 

C2 

O 

1 

U 

y 

= 

r;2 

r*2 

c3 

V 

z 

r3 

C1 

r3 

C2 

i 

Of) 

U 

w 

j 

(45) 


where 


pi  = 

sin  0 1 sin  a + cos  0 j sin  6w  cos  a 

1! 

*— CN 

U 

— sin  0 j cos  a + cos  0 | sin  5w  sin  a 

o 

ii 

cos  0^  cos  5w 

o 

h-N) 

II 

cos  5 w cos  a 

q2  - 

cos  5w  sin  a 

n 

U>N) 

II 

- sin  5w 

o 

II 

—cos  0^  sin  a + sin  0j  sin  6w  cos  a 

(-3  — 
c2 

cos  0 j cos  a + sin  0 j sin  5w  sin  a 

= 

C3 

sin  0j  cos  6w 

(46) 
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Figure  16.  Definition  of  the  Coordinate  Systems. 


TRACKER 


(b)  Rotation  About  the  V Axis. 


a - Angle  between  true  north  and  the  v axis  (radar 
coordinates)  measured  in  the  u,  v plane  posi- 
tive clockwise. 

6w  - Deflection  angle  of  the  local  gravity  vector 
at  the  radar  out  of  the  meridian  plane  posi- 
tive west. 

6u  - Deflection  angle  of  the  local  gravity  vector 
at  the  radar  in  the  meridian  plane  positive 
when  north  of  the  normal  to  the  reference 
ellipsoid. 

0*  - Geodetic  latitude  of  the  tracker,  positive  north. 

h - Height  of  the  tracker  above  the  geoid  (feet). 

N - Geoidal  separation-positive  when  the  geoid  is 
above  the  reference  ellipsoid. 

Figure  17.  Sequence  of  Rotations  Used  to  Derive 

the  Position  Transformation  Equations. 
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Each  rotation  matrix  in  Equation  (44)  is  orthogonal,  and  for 
this  reason,  the  inverse  of  Equation  (45)  is  simply: 


u 

q} 

X 

V 

= 

r*2 

r3 

c2 

y 

w 

_ 3 

c3 

c3 

°3_ 

z 

A"1  "=  AT 


(47) 


The  transformation  of  position  from  x,  y and  z to  earth 
centered  fixed  coordinates  £,  n and  c,  is  then  simply  the 
following  translation  (see  Figures  16  and  17a): 


i = (x  + xe) 

v = (y) 

f (z  + ze) 


(48) 


where 


xe  = rg  cos  6 + (h  + N)cos0  + 

z = r sin  0 + (h  + N)sin0* 
c c ' ' 


(49) 


rg  is  the  radius  of  the  reference  ellipsoid  at  the  latitude 
of  the  tracker,  which  is  given  by  the  following  equation: 


= a ( 1-e) 

V 1 - (2c  - c^)  cos^  0 


(50) 
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The  relation  between  geodetic  latitude  0*  and  geocentric 
latitude  0 may  be  easily  found  from  Equation  (49).  This  is 
as  follows: 


tan  0 - (1  - e)2  tan  0* 

Writing  Equation  (48)  in  more  convenient  form,  we  get 
f = x + cj2 

tj  = y 

? = z + c}2 

C|2  = re  cos  0 + (h  + N)  cos  0* 

C*2  = re  sin  0 + (h+N)sin0* 

a (1  -e) 

r — • 

^ 1 - (2e  - e2)  cos2  0 
0 = tan'l  f"(l  - e)2  tan  0*1 


(51) 


(52) 


(53) 


We  will  use  values  of  a and  e as  determined  by  the  Fisher 
Ellipsoid  (1960) 

a = 20,925,741  ft 

e = 3.35233051  x 10-3 


e is  ellipticity 
4. 5. 4. 3 Error  Propagation 

Using  the  equation  for  propagation  of  variance 

SX  - mAm*t 
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and  operating  on  Equation  (45)  to  find  the  matrix  of  partials 

M one  finds  that  M equals  the  matrix  of  C s. 

* x n 


Then  the  covariance  matrix  for  the  x,  y and  z coordinates  is 

S2  = ASjA7  (55) 

where  Sj  is  the  covariance  matrix  in  the  radar  pad  coordinates 
U,  V and  W. 

4.6  WHITE  SANDS  MISSILE  RANGE 


4.6.1  The  QD  Digital  Filter 

4.6. 1.1  In  early  1965,  due  to  expanding  requirements  for 
the  range  real-time  digital  computer,  the  need  for  a fast, 
efficient,  and  storage-conservant  data  filter  became  apparent. 
Mr.  W.  A.  McCool , Director  of  the  Analysis  and  Computation 
Directorate,  having  an  extensive  background  in  electronic 
engineering  and  analog  computing  developed  an  idea  and 
procedure  whereby  techniques  of  a basic  analog  feedback 
filter  could  be  applied  to  digital  data.  The  development  of 
concepts  of  real-time  data  filtering  at  WSMR  began  with  this 
original  idea,  and  the  credit,  although  others  have  contributed 
significantly  to  usage  development,  belongs  to  Mr.  McCool. 


The  first  operational  filter  developed  at  WSMR  under 
this  theory  was  the  Digital  Filter  Experimental  (DFX), 
and  it  was  used  operationally  as  late  as  December  1969.  DFX 
was  a fast  (less  than  1 millisecond  of  7044  time  per  component 
at  a 9-point  span)  and  efficient  (performed  equally  as  well 
as  constrained  least  squares  in  terms  of  smoothness  and  was 
more  conservant  of  computing  time  at  equivalent  point  spans) 
data  filter  used  for  all  types  of  range  safety  applications. 
DFX  was  employed  as  a variable  span  filter  with  changes  in 
span  occurring  according  to  observed  acceleration  in  the 
filtered  data.  The  computing  time  required  for  DFX  was 
directly  proportional  to  the  length  of  the  filter  span  and 
there  were  certain  elements,  in  array  form,  which  required 
retention  in  computer  core  storage. 

Many  ideas  innovated  in  DFX  carried  through  in  the 
enhancement  of  QD.  The  conversion  to  QD  was  implemented 
because  of  the  following  advantages:  (1)  filter  computing 
time  for  QD  is  the  same  regardless  of  span  (and  faster  than 
DFX  at  a minimum  equivalent  span);  (2)  core  storage  require- 
ments are  about  2/3  that  of  DFX;  and  (3)  QD  performance  over 
data  containing  random  noise  is  decidely  better  than  that 
of  DFX. 


The  QD  digital  filter  was  developed  at  WSMR  in 
1967,  and  it  has  since  been  used  in  computer  programs  for 
real-time  mission  support  and  in  a variety  of  other  types  of 
deferred  programs.  The  basic  concept  of  the  QD  filter  is 
not  new.  For  example,  the  alpha-beta  (a-8)  filter  was 
developed  more  than  ten  years  ago  for  target  position  predic- 
tion in  track-while-scan  radars.  The  significant  aspects  of 
the  QD  filter  lie  in  (1)  minimal  computing  time,  (2)  unique 
mathematical  formulation  based  on  an  nth  order  polynomial 
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curve  fit  and  (3)  extension  to  digital  data  smoothing.  The 
purpose  of  the  next  section  is  to  summarize  the  development 
of  the  QD  formulas  for  the  second  and  third  order  fits. 

4 . 6 . 1 . 2 QD  Theory 

Given  a set  of  discrete  data  values  with  equal 
increments  of  the  argument, 


x(tn)  = 


VX1 


2 ’ 


N ’ 


(1) 


h 


n+1 


(2) 


the  second  order  QD  formulas  for  sequential  filtering  these 
data  are  described  as  follows: 


QD  is  a recursive  predict-correct  procedure.  The  predic- 
tion formulas  are: 


x ..  = x , 
n+1  n ’ 


(3) 


where : 


x - second  derivative  estimate  from  preceding  step, 
n 


xr+1  - predicted  second  derivative; 


n+1 


xn  + hV 


(4) 


xn  + -^  - predicted  first  derivative, 


x - first  derivative  estimate  from  preceding  step; 
n 


— A /s  h 2 " 

x . . = x„  + hx  + q— x , 
n+1  n n 2 n 


(5) 
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or 


n+1 


= x_ 


2' xn+l 


xn> 


where 


(6) 


xR+1  - predicted  data  value, 

xn  - data  value  estimate  from  preceding  step. 


The  prediction  formulas  (4)  and  (5)  are  Taylor  series 
truncated  for  the  second  degree  curve  fit.  Formula  (6)  is  a 
more  computationally  efficient  version  of  formula  (4)  obtained 
by  substituting  formula  (4)  into  formula  (5).  Note  that 
formula  (6)  is  the  familiar  trapezoidal  integration  formula. 

e corrections  of  the  predicted  values  x , ,,  x .. 

n+1  n+1 

x are  based  upon  the  difference  between  the  given  and 

n-ri  D 

predicted  data  values,  i.e., 


Ax  = xn+l  - xn+l 
The  correction  formulas  are: 

A 

xn+l  = xn+l  +klAx 

A 

xn+l  =In+l  +k2Ax 
xn+1  = xn+l  + k3Ax 


(7) 


(8) 

(9) 

UO) 


where 


n+1 


"n+l 


current  estimate  of  the  true  data  value 

current  estimate  of  the  corresponding  first 
derivative 
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x , - current  estimate  of  the  correspond: ng 
second  derivative 

kj.kg.kg  - correction  coefficients  which  minimize 
the  error  of  the  estimates 

Formulas  (4),  (6),  (7),  (8),  (9)  and  (10)  prescribe  the 
computational  procedure  of  the  second  order  QD  digital 
filter.  The  QD  formulas  use  the  estimates  evaluated  in  the 
preceding  step  and  are,  therefore,  said  to  be  recursive. 

This  being  the  case,  QD  must  be  initialized,  i.e.,  initial 
estimates  must  be  obtained  independently.  The  most  common 
procedure  is  to  compute  the  initial  estimates  from  the  first 
few  data  values  with  a conventional  least-squares  curve  fit. 

The  QD  formulas  were  derived  from  the  constrained 
least-squares  (CLS)  filter  applied  to  M given  data  values 
with  estimates  determined  for  the  current  given  data  value, 
i.e.,  m=M.  In  the  CLS  filter,  intercept  and  slope  constraints 
are  applied  at  the  oldest  value  end  of  the  M value  span, 
i.e.,  m=0.  This  means  that  the  polynomial  fit  to  the  M data 
values  in  the  least-squares  sense  must  also  contain  the 
estimates  of  the  true  value  and  the  corresponding  first 
derivative  at  the  span  point  m=0.  Using  these  constraints 
in  the  second  order  CLS  filter,  the  predicted  data  value  and 
its  derivative  can  be  obtained  with  the  truncated  Taylor 
series . 


*n+l  ~xn-M+1  + Mh*n-M+1’ 

*n+l  = Xn-M+1  + M h VM+1  + ^7^“  Xn-M+1 


(11 

(12 


Where  the  subscript  n-M+1  corresponds  to  m=0  in  the  filter 

A 

span  and  ft  and  x are  the  constraints  which  were 

n-M+1  n-M+1 

computed  in  the  preceding  step  as  &n_M+2  anc*  *n-M+2  us*n6 


Taylor  series  again.  Since  S ^+j  and  ^n_^j+j  are  fixed, 

only  x can  be  adjusted  to  make  formula  (12)  fit  the  M 

n-M+1 

data  values  in  the  least-squares  sense.  Let  Ax  be  the 
required  correction.  Then 


x i»  •»- 1 Xn-M+1 


+ Ax. 


We  note  here  that  the  second  derivative  is  constant 
across  the  second  order  CLS  filter  span.  Applying  the 
correction  to  formulas  (11)  and  (12)  gives  the  estimates, 


A A A 

xti+l  = Xn-M+1  + M h (Xn-M+1  + Ax) 


A A 

xn+l  = x 


n-M+1  + Mh  x 


,(Mh)2 

‘n-M+1 


(5c*  + Ax). 


Subtracting  formulas  (11)  and  (12)  from  formulas  (14)  and 
(15),  respectively,  gives 


xn+l  = xn+l  + M h Ax- 


xn+l  " xn+l 


(Mh)~ 


Ax 


Comparison  of  formulas  (13),  (16)  and  (17)  with  formulas 
(8),  (9)  and  (10)  reveals  the  functional  relationship  among 
the  correction  coefficients  k^,  kg  an£i  kg.  From  formulas 
(13)  and  (8)  we  have 

Ax  = k j Ax. 


Then  from  formulas  (16),  (9)  and  (18), 
k2  = M h kr 


— pwpp  ’ 


(13) 


(14) 

(15) 

(16) 

(17) 

(18) 

(19) 


98 


and  from  formulas  (17),  (10)  and  (18), 


k3  = iMhl2  k 


(20) 


In  the  QD  theory  k1  is  developed  as  a function  of  the  corres- 
ponding CLS  filter  span  M,  i.e., 

kl  = (21 ) 

Thus,  when  M is  arbitrarily  specified,  the  QD  correction 
coefficients  are  determined  with  formulas  (21),  (19)  and 
(20). 

We  note  that  QD  prediction  formulas  4 ) and  (5)  do  not 
contain  the  constraints  $n_M+1  and  *n_M+1  which  are  eliminated 
in  the  QD  theory.  This  is  easily  shown  as  follows.  Recog- 
nizing that  any  two  points  contained  in  the  fitted  curve  of 
the  preceding  step  can  be  related  by  a Taylor  series,  we  can 
write 


A A A 

xn  = VM+1  +(M-l)h  xn.M+1 
xn  = VM+1  + Vm+I 

Xn-M+1 


(22) 

(23) 


QD  prediction  formulas  (4)  and  (5)  are  obtained  by  combining 

formulas  (22),  (23),  (11)  and  (12)  to  eliminate  & and 

* n-M+l 

*n-M+l'  ^us  ■ * t i-s  now  clear  that  the  second  order  fitted 
curve  established  in  the  preceding  step  and  containing  the 
intercept  and  slope  constraints  is  used  for  the  prediction 
phase  of  the  current  computing  step,  i.e.,  formulas  (3),  (4) 
and  (6). 
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The  computing  formulas  of  the  third  order  QD  digital 
filter  are  as  follows: 


xn+l  = xn  + h*n 

(24) 

*n+1  = xn  +y-  (Vn+*n) 

(25) 

x ='x'  +.!l_  (T  + x ) 

n+1  n 2 'ii  + l iv  1 2 

(26) 

Ax  = xn+1  -xn+l 

(27) 

*n+1  =\  + klAx 

(28) 

A /V 

xn+l  = xn  + k2Ax 

(29) 

A A 

xn+l  = *„  + k3Ax 

(30) 

xn+l  = x„  + k4Ax 

(31) 

In  the  third  order  QD  formulation  the  intercept  and 
slope  constraints  are  applied  and  the  second  and  third 
derivatives  are  adjusted  to  satisfy  the  least-squares  error 
criterion.  The  correction  coefficients  k^,  k?  , k^  and  k^ 
are  developed  as  functions  of  M,  the  equivalent  third  order 
CLS  filter  span.  Note  that  formulas  (25)  and  (26)  have  been 
modified  for  computing  efficiency. 

In  its  basic  formulation,  QD  is  a real-time  filter, 
i.e.,  the  argument  of  the  estimates  corresponds  to  the 
latest  value  accepted  by  the  filter.  On  the  other  hand, 
estimates  can  be  obtained  for  the  data  value  at  the  oldest 
end  of  the  span  which  are  significantly  better  than  the 
corresponding  real-time  estimates.  Such  estimates,  which 
are  often  called  "smoothed"  data,  correspond  to  the  constraints 
computed  in  the  CLS  filter.  Smoothed  estimates  are  simply 
obtained  with  a Taylor  series  expanded  about  the  real-time 
estimates , i.e.. 
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for  the  third  order  QD  filter. 
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The  basic  QD  filter,  used  for  WSMR  safety  applications, 
operates  at  variable  point  spans  which  are  adjusted  according 
to  acceleration.  The  filter  spans  vary  from  as  few  as  17  to 
as  many  as  46  points.  When  QD  is  used  to  filter  FPS-16 
data,  2 raw  data  points  are  averaged  to  form  1 filter  input 
point.  The  span  of  the  filter  is  varied  according  to  a 
philosophy  that  the  higher  the  acceleration  the  shorter  the 
span.  Thus,  when  accelerat  io..  j are  high,  maximum  filter 
response  to  changes  in  the  input  curve  are  achieved;  when 
accelerations  are  low,  the  benefits  of  maximum  smoothing  are 
available. 


The  filter  is  always  initialized  or  re-initialized  at 
the  minimum  span.  Should  span  adjustment  be  deemed  necessary, 
the  span  is  lengthened  or  shortened  by  only  1 point  per 
filter  cycle.  A table  of  acceleration  values  determines  what 
span  the  filter  should  be  or  whether  the  span  should  be 
increasing  or  decreasing.  This  table  of  values  can  be  formu- 
lated according  to  the  performance  of  a particular  vehicle 
to  both  minimize  "overshoot"  and  achieve  maximum  smoothing 
as  necessary . 
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4.6. 1.3  Outstanding  Features 


Several  features  presented  here  will  make  WSMR  QD 
filter  operational  use  more  clear  and  will  explain  why  it  is 
extensively  applied. 

a.  QD  Spike  Filter.  QD  has  an  inherent  spike 
filter,  not  previously  discussed,  which  vastly  improves 
filter  operation.  The  recursive  nature  of  the  filter  is 
employed  such  that  after  the  prediction  equation  is  applied 
and  AX  is  computed  (AX  = input  position  minus  predicted 
input  position),  the  magnitude  of  AX  is  compared  to  a pre- 
determined spike  limit  L.  L is  an  amount  of  position  change 
reflecting  an  acceleration  that  the  particular  vehicle  under 
observation  is  not  capable  of  achieving.  L must  be  formulated 
such  that  values  of  AX  equal  to  or  exceeding  L can  be 
considered  to  be  noise  or  spikes  in  the  raw  data. 

When  a value  of  AX  is  observed  to  equal  or 
exceed  L,  a spike  counter  (C)  is  incremented  and  the  predicted 
value  of  position  (Xn  + 1)  is  used  for  filter  input  rather 
than  the  raw  input  position  (Xn  + 1).  This  substitution  is 
allowed  to  occur  n/2  times  for  a given  filter  span  (n).  If 
n/2  + 1 spikes  are  observed  on  a filter  cycle,  re-initializa- 
tion is  automatically  forced.  In  general  usage  at  WSMR,  a 
double  spike  limit  is  employed.  If  | AX | > LI  and  Cl^n/2  + 1, 
then  the  same  test  is  made  for  an  L2  and  C2 . The  Cl  and  C2 
counters  are  incremented  and  decremented  according  to  the 
magnitude  of  AX.  For  example,  if  the  value  of  AX>L1  but 
<L2,  only  Cl  is  incremented.  This  double  spike  limit  allows 
a 'ine  degree  of  spike  editing  using  LI  with  L2  available  to 
delete  larger  spikes  and  avoid  unnecessary  re-initialization. 

b.  Initialization.  A 5 point  least-squares 
curve  fit  is  used  for  initialization  of  the  recursive  QD 
filter;  the  5 point  initialization  array  is  shifted,  and  the 
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latest  raw  data  value  is  inserted  on  each  cycle  of  the 
filter.  This  can  enable  the  filter  to  be  initialized  in 
less  than  5At  (At  being  the  data  sampling  interval).  In  the 
case  of  meteorological  use  of  the  QD  filter  for  computation 
of  wind  velocity,  a 7 point  least-squares  curve  fit  is  used 
for  both  initialization  and  to  form  the  input  to  the  filter 
on  each  cycle.  This  is  opposed  to  using  the  raw  input  data 
for  filter  input  and,  when  used  with  the  midpoint  smoothing 
technique  explained  in  subparagraph  4.6. 1.2,  provides  optimum 
smooth  data  which  lags  real  time  by  approximately  2 seconds. 
Studies  are  currently  underway  to  modify  the  meteorological 
QD  for  use  as  a real-time  filter  for  filtering  rough  singular 
measurements.  This  would  be  advantageous  for  forming  multiple 
station  solutions  from  some  instrumentation  previously 
described . 

c.  Variable  Filter  Span.  Variable  filter  spans 
allow  more  general  use  of  the  filter  through  the  increased 
versatility  gained  by  automatically  adapting  the  filter  to 
changes  in  observed  acceleration. 

d.  Filter  Speed  and  Storage  Requirements.  QD  is 
a real-time  data  filter  that  precisely  meets  the  design  and 
development  goals.  The  minimum  storage  requirements  and 
computing  speed  of  QD  must  be  rated  as  outstanding  features. 

4. 6. 1.4  Limitations 

Since  other  ranges  would  undoubtedly  have  to 
adjust  the  QD  filter  to  satisfy  their  particular  needs,  such 
adjustments  could  be  classed  as  limitations. 

The  only  limitation  or  possible  problem  area  at 
WSMR  with  QD  is  that  extremely  rough  input  data,  which 
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passes  through  the  spike  filter,  will  be  interpreted  as  an 
acceleration  change  and  will  cause  undesired  shortening  of 
the  filter  span.  While  this  is  not  a major  problem  at  WSMR, 
manual  overrides  have  been  programmed  into  the  WSMR  QD 
filter  for  fixing  the  filter  at  a particular  span  length 
and  for  providing  different  values  of  the  spike  limits. 

This  was  done  primarily  for  orbital  support  operations 
and  aircraft  tracking,  rather  than  range  safety  during 
missile  flight  experiments. 

4.7  ARMAMENT  DEVELOPMENT  AND  TEST  CENTER 

4.7.1  Least-Squares  Polynomial  Moving  Arc  Filter  Using 

Recursive  Sums 

The  real-time  range  safety  filter  in  current  use 
at  the  ADTC  employs  a modified  recursive  scheme  for  computing 
the  sums  required  to  perform  least-squares  polynomial  moving 
arc  smoothing. 

The  matrix  equation  for  this  type  of  smoothing  may 
be  written  as: 


A = C 1B 


where  A is  the  vector  of  polynomial  coefficients,  C is  the 
normal  matrix  with  elements 


n 

E 

m=l 
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i+j-2 


i = 1,2 d+1 

j = 1,2 d+1 


and  B is  a vector  whose  kth  element  is  of  the  form 
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n 

5:  y 

t=i 


for  k = 0,1,2 d, 
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d being  the  degree  of  the  polynomial;  t^,  the  time  of  the 
sample  referenced  to  the  midpoint  of  the  span  under  considera- 
tion; and  y^,  the  coordinate  (x,y  or  z)  at  time  t^. 

For  instance,  if  one  has  obtained: 


biO)  - I y^r1'  - Wi 

k-1  2 


and  wishes  to  find: 
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one  may  do  so  in  terms  of  b.(1)  by  means  of  the  recursion 
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The  elements  of  C are  computed  similarly  by  considering  the 
yk=l.  However,  if  the  data  are  evenly  spaced,  C is  uniquely 
determined  when  the  degree  of  the  polynomial  and  the 
number  of  points  smoothing  has  been  selected,  and  C-1  can  be 
pre-computed  and  stored  for  use  in  the  smoothing  computations. 
In  this  case,  A=1  and  the  equations  for  recursion  are  simpli- 
fied. With  this  formulation  no  lengthy  summations  are  required, 
and  the  number  of  operations,  once  initiated,  are  independent 

of  the  number  of  points  used  in  the  smoothing.  Furthermore, 

o 

by  computing  recursively  E y^,  one  can  obtain  the  sum  square 
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residuals  of  each  span  as  an  estimate  of  the  error  in  the 
data  by 


S = T.  yj  - ATB 
•'k 

4.7. 1.1  Features.  The  ADTC  filter  algorithm  is  recursive 
in  the  sense  that  it  uses  the  output  of  a previous  fixed 
time  span  as  input  to  the  next  time  span.  Data  at  the 
newest  time  point  is  used  to  modify  the  preceding  estimate, 
and  the  contribution  of  data  at  the  oldest  point  in  the 
preceding  span  is  deleted  from  the  estimate.  This  method, 
once  the  filter  is  initiated,  allows  for  rapid  calculation 
and  minimal  computer  core  requirements,  and  leads  to  a 
simple  estimate  of  individual  sensor  data  quality  for  use  in 
source  selection. 

4.7. 1.2  Limitations.  The  algorithm  is  unstable  in  the 
sense  that  it  is  unable  to  eliminate  the  effect  of  any  past 
erroneous  data  on  the  recursive  estimations.  Filter  initiali 
zation  can  be  lengthy.  For  initialization  of  a 51-point 
filter,  at  a data  sample  rate  of  10  samples  per  second.  5 
seconds  is  required.  The  filter  is  most  efficient  for  track 
of  objects  with  smooth  trajectories  and  is  not  adaptive  to 
those  performing  intricate1  maneuvers. 

4.7.2  Planned  Filter  Applications 

The  existing  ADTC  filter  moot,  was  initially 
developed  for  a post-flight  data  reduction  application. 

When  used  in  this  way,  lengthy  initialization  times  and  lack 
of  adaptiveness  and  stability  were  not  of  serious  consequence 
The  adaptation  of  the  filter  to  real  time  did  not  overcome 
these  problems.  As  a result,  the  present  filter  is  most 
effective  on  smooth  flight  profiles  with  radar  data  of  good 
tracking  quality.  However,  it  is  common  for  the  tactical 
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type  of  vehicles  tracked  by  ADTC  radars  to  perform  intricate 
maneuvers  at  low  elevation  tracking  angles.  The  tracking  of 
more  than  one  vehicle  during  a mission  is  also  common  as  is 
the  tracking  ot  rapidly  accelerating  air-to-air  and  air-to- 
ground  missiles  launched  from  maneuvering  aircraft.  The 
ADTC  is  presently  in  the  midst  of  developing  a centralized 
CRT  oriented  control  center  capability  to  replace  existing 
decentralized  plotboard  facilities.  Effective  use  of  the 
new  concepts  and  technology  require  that  the  existing  computer 
software,  algorithms  and  procedures  be  totally  revamped.  To 
adequately  support  these  situations,  a filter  for  specific 
ADTC  application  has  been  designed.  (References  15  and  16.) 

4.7.2. 1 Filter  Method.  The  replacement  for  the  current 
ADTC  real-time  filter  utilized  a multiple  exponential  filter 
technique  for  each  data  channel.  The  design  application  has 
the  following  features: 

a.  Accepts  data  from  multiple  radars. 

b.  Filter  weights  may  be  selected  to  give  zero,  first 
or  second  degree  polynomial  filtering  with  QD,  Kalman  or  no 
constraints  imposed. 

c.  Three  sets  of  filter  weights  may  be  applied  to 
each  channel  of  input. 

d.  A zero  degree  polynomial  may  be  used  to  estimate 
noise  variance. 

e.  Adaptation  is  accomplished  through  inter-filter 
velocity  comparisons  and  residual  trend  detection. 

f.  Initialization  is  accomplished  by  means  of  an 
expanding  memory  least-squares  polynomial  algorithm. 


The  filter  scheme  makes  use  of  the  fact  that  several  expon- 
ential filters  can  be  applied  in  less  computer  time  than  is 
required  for  the  present  filter.  Each  data  channel,  then, 
is  assigned  three  different  exponential  filters  representing 
different  effective  memories  or  bandwidths.  The  calculations 
are  performed  for  each  filter  at  each  data  point.  The 
optimum  filter  for  the  prevailing  conditions  of  data  quality, 
vehicle  maneuvering  and  radar  performance  is  selected  through 
inter-filter  comparisons  for  that  channel.  This  enables  the 
filters  to  adapt  to  prevailing  conditions  to  avoid  long 
settling  periods  after  abrupt  maneuvers,  missile  separations 
and  data  dropouts. 
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